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ABSTRACT 

Using the public data from the Herschel wide field surveys, we study the far-infrared properties of 
optical-selected quasars from the Sloan Digital Sky Survey. Within the common area of ^ 172 deg^, 
we have identified the far-infrared counterparts for 354 quasars, among which 134 are highly secure 
detections in the Herschel 250 pm band (signal-to-noise ratios > 5). This sample is the largest far- 
infrared quasar sample of its kind, and spans a wide redshift range of 0.14<z < 4.7. Their far-infrared 
spectral energy distributions, which are due to the cold dust components within the host galaxies, 
are consistent with being heated by active star formation. In most cases (> 80%), their total infrared 

luminosities as inferred from only their far-infrared emissions already exceed 10^^ Lq, and thus 

these objects qualify as ultra-luminous infrared galaxies. There is no correlation between and 

the absolute magnitudes, the black hole masses or the X-ray luminosities of the quasars, which further 
support that their far-infrared emissions are not due to their active galactic nuclei. A large fraction of 
these objects (> 50-60%) have star formation rates > 3OOM0yr“^. Such extreme starbursts among 
optical quasars, however, is only a few per cent. This fraction varies with redshift, and peaks at around 
z ^ 2. Among the entire sample, 136 objects have secure estimates of their cold-dust temperatures 

(T), and we find that there is a dramatic increasing trend of T with increasing We interpret 

this trend as the envelope of the general distribution of infrared galaxies on the (T, L^j^) plane. 

Subject headings: infrared: galaxies; galaxies: starburst; galaxies: high-redshift, galaxies: evolution, 
(galaxies:) quasars: general 


1. INTRODUCTION 

Ultra-Luminous InfraRed Galaxies (ULI RGs), discov¬ 
ered as a distinct class in the early 1980 ’s (jHonck et al.l 
11984 Il985l : lAaronson fc 01szewskilll98l ) by the Infrared 
Astronomical Satellite (IRAS) survey of the sky in 12 
to 100 pm, are believed to host extreme star formation 
regions that are heavily enshrouded by dust. They are 
characterized by their exceptionally high IR luminosi¬ 
ties {Ljr > 10^^ Lq; integrated over restframe 8 to 
1000 pm), which are believed to be predominantly due 
to the re-radiation o f star light processed by dust (see 
iLonsdale et al.l [20061 for a review), and thus imply very 
high star formation rat es (SFR) of > 10 0-1000 
(using the conversion of lKennicuttlll998h completely hid¬ 
den by dust. 

However, it has been noticed ever since their discovery 
that a significant fraction of ULIRGs, especially those 
very luminous ones, have optical signat ures ind i cative 
of classic AGN activities. For examples, iCarterl (|1984f ) 
notes that ten among a sample of 13 IRAS sources with 
60 pm flux density /bo > 1.2 Jy are Seyfert galaxies. 
iSanders et al.l (|1988[ ) show that ten galaxies among the 
324 sources with /eo > 5.4 Jy from the IRAS Bright 
Galaxy Survey are ULIRGs and that they all have a 
mixture of starburst and AGN signatures, which has led 
them to propose an evolutionary scenario that ULIRGs 
are the prelude to quasars. In their redshift survey of 
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the IR AS >560 > 0.6 Jy galaxy sample, iLawrence et al.l 
(|1999f ) have found that ^ 20% of their 95 ULIRGs are 
AGN. The consensus is that ULIRGs that have “warm” 
IR colors, i.e., whose emissions tend to peak at restframe 
mid-IR (MIR) rather than far-IR (FIR), generally host 
(optical) AGN, which could be the main energy sources 
that p ower their strong IR emissions (e. g . Jde Griip et al. 
1985 : lOsterbrock fc De RobertisI II985I: iKim fc Sanders 
1998). On the other hand, “cold” ULIRGs that have 
their IR emissions pe ak at FIR should be mostly powered 
by st arbursts (e.g., lElston et al.l 119851: iHeckman et al.l 

[iM3). 

An interesting question then is whether any AGN 
ULIRGs, especially quasar ULIRGs, have starbursts that 
dominate their IR emissions. Quasars represent the most 
extreme process of supermassive black hole accretion, 
while starbursts are the most extreme process of star for¬ 
mation. It is expected that the interplay of these two ex¬ 
tremes will have important consequences. This is made 
particul arly important by th e quasar evolutionary sce¬ 
nario of ISanders et al.l (| 19881) . as such objects could be 
the transitional type between non-quasar ULIRGs and 
“fully exposed” quasars. Sanders et al. themselves be¬ 
lieve that AGN heating is the main mechanism for the 
strong IR emission of ULIRG quasars. In their dis¬ 
cussion of PG quasar continuum di s tribut ions from UV 
to millimeter (mm), ISanders et al.l (| 19891) further pro¬ 
pose that a warped galactic disk (beyond the central 
^ 10 pc to a few kpc) heated by the central AGN can 
explain the ent ire range of IR emissio n from 5 pm to 
1 mm. However. iRowan-RobinsonI (|1995f ) argues that this 
scheme could only be viable when the total luminosity 
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in IR is comparable or less than that in UV-optical; in¬ 
stead, he has successfully modeled the IRAS-detected PG 
quasars by attributing their mid-IR emissions to AGN 
heating and the FIR emission to starburst heating, re¬ 
spectively. More analysis using larger samples from the 
Infrared Space Observa tory (ISO) observa tions support 
this view. For example, iHaas et ^ (j2QQ3f ) have studied 
64 PG quasars and have concluded that starburst heat¬ 
ing is more likely the cause of the observed cold dust 
30-50 K) FIR emissions among the ULIRG part of 
the sample. However, they have also pointed out that 
AGN heating should be the main power source for those 
extreme ones that qualify as “Hyper-luminous Infrared 
Galaxies” (HyLIRG; usually defined by Ljr > 10^^ Lq). 

There are also two pieces of important, albeit in¬ 
direct, evidence supporting that the FIR emissions of 
quasar ULIRGs are likely due to starbursts. First, a 
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a strong indicator of active star formation. Second, 
most quasar ULIRGs have polycyclic aromatic hydro¬ 
carbon (PAH) features, which are also strongly indica¬ 
tive of on-going star formation (e.g.. iSchweitzer et ail 
aJM et al. I 2 OO 7 I: iHao et aT]l200'ii: iNetzer et al.ll2007[ 


Gao et aT]l2008 b While none of these are sufficient to as¬ 


sert that starbursts dominate the strong FIR continua of 
quasar ULIRGs, it is clear that they at least contribute 
significantly. 

The picture above is largely based on the ULIRGs in 
the nearby universe where they can be studied in de¬ 
tail. It would not be surprising if any of it changes 
at high redshifts, as both quasars and ULIRGs evolve 
strongly. The number density of quasars rises rapidly 
from z = 0 to z = 1, and reaches the peak at z ^ 2-3 
(e.g.. lOsmerl 120041 ). Similarly, while ULIRGs are very 
rare objects today, they are much more numerous in 
earlier epochs. The deep ISO surveys have revealed 
a large number of IR-luminous galaxies, among which 
> 10% are ULIRGs and m any are at 2 : > 1 (e.g., 
lR.owan-R.obinson et al.l[2QQ3 ). The discovery of the so- 
called “snbmillin ieter (snbmm) ga laxies” (SMGs) at 450 
and 850 pm (see iBlain et~al]l2002l for a review) added a 
new population to the ULIRG family, as most of them 
are at z ~ 2-3 and have Ljr > 10^^ Lq dominated by 
the emission from cold dust. Furthermore, Spitzer ob¬ 
servations suggest that high stellar mass (> 10^^ Mq) 
and otherwise “normal” st ar-forming galaxie s at z ^ 2 
are likely all ULIRGs (e.g.. lDaddi et alj2005f ). which in¬ 
creases the ULIRG number density at high redshifts to a 
more dramatic level than expected. On the other hand, 
molecular gas has also been detected in qnasy UL IRGs 
from 2 : > 1 to 6 (jSolomon fc Vanden BontI 12005 . and 


the re ferences therein; see also e.g., I Wang et ahl l2010l 
l2oTT^ for the recent results at 2 : ^ 6), lending sup¬ 
port to the starburst-powered interpretation of the FIR 
emission of such objects at high redshifts as well. 

If quasars and dust-enshrouded starbursts do co-exist, 
it is important to investigate their co-evolution, which 
would require a large sample over a wide redshif t range. 
Herschel Space Observatory (jPilbratt et al.l[20Tol ) has of¬ 
fered an unprecedented opportunity to investigate this 
problem at the FIR wavelengths that were largely unex¬ 
plored by the previous studies. Herschel had two imag¬ 


ing instruments, namely, the Photo detector Array G am- 
era and Spectrometer (PAGS, iPoglitsch et al.l[2QlQf ) and 
the Spectral and P hotometric Imaging REceiver (SPIRE, 
IGriffin et, a,l .1 l2(rfnh . The PACS bands are 100 (or 70) 
and 160 pm, and the SPIRE bands are 250, 350 and 
500 pm. Together they sample the peak of heated dust 
emission from 2 ) ~ 0 to 6 and beyond. There already have 
been a number of studies on the EI R emission of quasars 
using Herschel observati o ns (e.g.. iSerieant et al.l 120101: 
Leipski et al.l l2010L 120131: iDai et al.l 120121: iNetzer et al l 
20141 ). however the current collection of quasars that 
have individual Herschel detections are still very scarce 
in number and few have spanned a su fficient redshift 
range (for exam ples. iLeipski et al.l (j2013f ) present 11 ob¬ 
jects at 2 ) > 5: iDai et al.l (l2012D i nclude 32 objects at 
0.5 < 2 ) < 3.6; INetzer et al.l (j201^ report ten within a 
narrow window at 2 : ~ 4.8). 

In this paper, we present a large sample of optical 
quasars that are detected by the Herschel^ and provide 
our initial analysis of their EIR properties. The quasars 
are fr om the Sloan Digital Sky Survey (SDSS: lYork et aP 
[ 2 OOOI) . and the EIR data are from the public releases 
of four major wide field surveys by Herschel^ namely, 
the Herschel Astrophysical Tp ahertz Large Area Survey 
(H-ATLAS: lEales et al.l[20Tol ). t he Herschel Multi-t iered 
Extragalactic Survey fHer MES: lOliver et al.l[20T^ . the 
Herschel Stripe 82 Survey (jViero et al.l 2014 HerS;), an d 
the PAGS Evolutionary Probe (PEP: iLntz et al.ll20Ilf ). 
We describe the data and the sample construction in ^ 
and present our analysis of the EIR dust emission in ^ 
The implications of our results are detailed in 21 we 
conclude with a summary in 21 The catalog of our sam¬ 
ple is available as online data in its entirety. All quoted 
magnitudes in the paper are in the AB system. We 
adopt the following cosmological parameters throughout: 
Qm = 0.27, Da = 0.73 and Hq = 71 kms“^ Mpc“^. 

2. DATA DESCRIPTION AND SAMPLE CONSTRUCTION 

In brief, we built our sample by searching for the coun¬ 
terparts of the SDSS quasars in the Herschel wide field 
survey data. Eor the sake of simplicity, hereafter we refer 
to these objects as “IR quasars”. We describe below the 
data used in our study and the constructed IR quasar 
sample. 


2.1. Parent quasar samples 

The parent quasar samples that we used are based on 
the SDSS Data Release 7 and 10 quasar catalogs (here¬ 
after DR7Q and DRIOQ, respectively), which are sum¬ 
marized as follows. 

DR7Q: As detailed in iSchneider et al.l (j2010f ). this 
quasar catalog is based on the SDSS DR7. It con¬ 
cludes the quasar survey in the SDSS-I and SDSS-H 
over 9380 deg^, and supersedes all previously released 
SDSS quasar catalogs. It includes 105 783 quasars be¬ 
tween 2 ) = 0.065 and 5.46 (the median at 2 : = 1.49), 
all with absolute i-band magnitudes (M^) brighter than 
—22 mag. 

DRIOQ: This quasar catalog is derived from the on¬ 
going Baryon Oscillation Spectroscopic Survey (BOSS) 
as part of the SDSS-HL While it was released in the 
SDSS DRIO, its main target selection was based on the 
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Figure 1. Comparison of the photometry, the redshifts and the absolute magnitudes from the DR7Q and the DRIOQ for the overlapped 
population in these two quasar catalogs (16 356 objects in total). For this particular subset, we adopt the DRIOQ values in this paper 
unless noted otherwise. For the sake of consisten cy, we use t he absolute magnitudes ^-corrected to ^ = 2 for the DR7Q quasars as well, 
adapting the values from the catalog presented in IStien et al ] (IMTlIb 



Figure 2. Redshift distribution of the SDSS DR7Q (red curves) 
and DRIOQ catalogs (yellow curves). The solid ones are for the 
entire catalogs, while the dashed ones are for the “ho mogeneous” 
subsa mples as described in ISchneider et ahl (|2010l j and IParis et al.l 

SDSS DR8. The detailed desc ription of the catalog can 
be found in IP aris et al.l (|2Q14D . It includes new quasars 
with Mi[z = 2] < —20.5mag from the SDSS-III, where 
Mi [z = 2] is the a bsolute magnitud e /^-corrected to z = 
2 (for details, see IP ar is et al.l [201^ . It also includes a 
large number of known quasars of similar characteristics 
(mostly from SDSS-I and II) that were re-observed by 
BOSS. In brief, the catalog contains 166 583 quasars 
over 6373deg^, with redshifts ranging from 0.05 to 5.86. 


data set over the widest area to date. The redshift 
distributions of DR7Q, DRIOQ and the merged cata¬ 
log are shown in Figure [2l For simplicity, hereafter we 
refer to the merged sample as the “SDSS quasar sam¬ 
ple” and the quasars therein as the “SDSS quasars”. We 
note that DR7Q and DRIOQ are statistically different 
samples, and therefore any statistical results from this 
merged catalog should be inferred with caution. This is 
exacerbated by the fact that the quasars were not se¬ 
lected u niformly within e i ther p R7 Q or DRIP 


Q or DRIOQ , as de - 
tailed in iSchneider et al] (l2(rfnh andlPlriseEaD (l2(rfl) . 
respectively. For example, while most of the quasars in 
D R7Q were selected usi ng the algorithms as described 
in iRichards et al.l (|2002f ) (with USELFLAG=1 in the cat¬ 
alog), a non-negligible fraction of them were selected 
early in the SDSS campaign when such algorithms had 
not yet been fully developed. DRIOQ is even more 
non-uniform in this sense, because only about half of 
its quasars (called the “CORE” sample) were selected 
uniformly (with U NIF0RM>0 in the c atalog) through the 

and the other half 


XDQSO method 
(called the “BONUS” sam ple) were selected using vari¬ 
ous different methods (see iRoss et al.l [20T^ . Neverthe¬ 
less, the non-uniformity does not affect our current work 
since we limit our study to the FIR properties of optical- 
selected quasars, whose being selected did not use any 
FIR information and thus should not favor or against 
any given FIR property. 


The quasars from these two catalogs are largely in¬ 
dependent, however there are still 16 356 of them being 
duplicates, which we define as the ones falling within a 
matching radius of 0.4". Figure [T] shows the compari¬ 
son of the photometry, the redshift measurements and 
the absolute magnitudes from these two catalogs for this 
overlapped population, which all agree reasonably well. 
For the sake of simplicity, we adopt the DRIOQ values 
in this work for these duplicates unless noted otherwise. 

In the end, we produced a merged sample of 256 010 
unique quasars, which represents the largest optical 
quasar sample selected based on the most homogeneous 


2.2. Herschel Data 

We utilized all high-level, publicly available data from 
the wide field Herschel surveys. In all cases we adopted 
the latest catalogs released by the survey teams to con¬ 
struct the FIR spectral energy distributions (SEDs). The 
basic characteristics of these surveys are summarized in 
Table [H and their relevant details are briefly described 
below. 

2.2.1. HerMES 

HerMES was the largest Herschel guaranteed time key 
program, and surveyed 100 deg^ in six levels of depth 
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Table 1 

Summary of Herschel wide field data and IR QSO sample 


Survey 

Field 

Level 

Goverage 

(deg 2 ) 

5 (7 limit 
(mJy) 

All^ 

Match 
< 3" 

SNR 
> 5 

*8250 > 56.6 
mJy 

HerMES 

COSMOS'^ 

LI 

5.00 

31.4 

216 

47 

7 

3 


GOODS-North ^ 

L2 

0.63 

27.1 

17 

4 

1 

1 


Bootes-HerMES 

L5 

11.29 

31.3 

426 

50 

18 

2 


EGS-HerMES ^ 

L5 

3.12 

29.7 

53 

6 

3 

0 


(Groth-Strip) 

L3 


30.9 


3 

2 

0 


ELAIS-NI-HerMES 

L5 

3.74 

30.4 

22 

3 

3 

0 


Lockman-SWIRE ^ 

L5 

19.73 

37.0 

276 

35 

12 

5 


(Lockman-East-ROSAT) 

L3 


31.5 


1 

0 

0 


(Lockman-North) 

L3 


31.4 


2 

1 

0 


ELS 

L 6 

7.30 

35.4 

89 

17 

5 

3 


XMM-LSS-SWIRE 

L 6 

21.61 

48.5 

402 

25 

3 

3 


(UDS) 

L4 


33.8 


3 

1 

0 


(VVDS) 

L4 


33.7 


8 

2 

0 

HerS 

Stripe82 

-L7 

80.76 

51.8 

4519 

132 

58 

53 

HATLAS 

SDP 

N/A 

19.28 

34.1 

690 

18 

18 

12 

Total 



172.46 


6710 

354 

134 

82 


Note. — (1) The detection limits are based on the SPIRE 250 pm source count histograms where the count drops to 50% of the peak 
value. The covered areas are calculated based on the pixel counts in the 250 pm maps. ( 2 ) The HerMES program has several “nested 
fields”, which are the sub-fields of deeper observations embedded in the shallower but wider parent fields. In this table, the nested fields 
are labeled by parenthesis and are placed after their parent fields. The numbers of matched quasars in these nested fields and the wider 
regions outside are quoted separately. 

^ “All” lists the number of SDSS quasars fall within the covered region. 

^ The PEP data are available in these fields. 


and spatial coverage combinations (“LI” to “L6”), using 
both SPIRE and PACS. Its latest data release, “DR2”, 
contains the SPIRE maps and the source catalogs in all 
three bands (|Wang et al.ll2Q13h . The PACS data have 
not yet been released. In this work, we used the wide- 
held data and excluded those galaxy cluster helds. 

The HerMES DR2 includes two sets of source catalogs 
based on different methods, one using the SUSSEXtrac- 
tor (SXT) point source extractor and the other using the 
iterative source detection tool StarEinder (SE) combined 
with the De-blended SPIRE photometry algorithm (DE- 
SPHOT). We adopted the band-merged version of the 
latter (denoted as “xID250”in DR2), which has the DE- 
SPHOT multi-band (250, 350 and 500 pm) photometry 
at the positions of the SE 250 pm sources. The 5 cr detec¬ 
tion limits in 250 pm range from ^ 30 to 48mJy in our 
helds of interest. 

2.2.2. H-ATLAS 

H-ATLAS was the largest Herschel open-time key pro¬ 
gram. It surveyed 550 deg^ using both SPIRE and PACS. 
Currently, this program has released the image maps 
and the source catalogs in the held observed during the 
Herschel Scienc e Demonstration Phase (H-ATLAS SDP ; 
llbar et ahl 120101: iPascale et al.|[2QTTI: iRigbv et al.lDOllT) . 
which covers ^19.3 deg^ and has reached the 5 a limit of 
33.4 mJy in 250 pm. 

The catalog that we adopted is the band-merged one 
with both the SPIRE and the PACS photometry. Eor 
the SPIRE bands, the source extraction was done using 
the Multi-band Algorithm for source extraction (MADX; 
Maddox et al. in prep.), which employed a localized back¬ 
ground removal and PSE hltering procedures to the map 


and extracted the sources using the 250 pm positions as 
the priors. Eor the PACS bands, aperture photometry 
was performed on the 100 and 160 pm maps at the SPIRE 
250 pm positions. 

2.2.3. HerS 

HerS observed r\j SO deg^ in the SPSS S tripe 82 region 
(|Abazaiian et al.1120091: lAnnis et al1l2014D using SPIRE, 
reaching the nominal 5 cr limit of 50.4 mJy in 250 pm. 
Both the image maps and th e source catalogs have been 
released (|Viero et all I2014D . We adopted the band- 
merged catalog, which was based on the SE 250 pm detec¬ 
tions and the DESPHOT photometry (|Roseboom et al.1 

[2oToh . 

2.2.4. PEP 

PEP was also a Herschel guaranteed time key program, 
which used PACS to survey six well-studied extragalactic 
fields and also a number ojf galaxy clusters. Both the im¬ 
age maps and the source catalogs have been made public 
through Data Release 1 (DRl) of the team. These data 
are not listed in Table [H All the PEP fields are within 
the HerMES fields, however they only covered a small 
fraction. 

Whenever possible, we used their 100 and 160 pm 
measurements in the COSMOS, GOODS-North, EGS, 
Lockman Hole fields to supplement the HerMES SPIRE 
photometry to better constraint the EIR SEDs of the 
detected quasars. These measurements were taken 
from the SE “blind e xtraction” catalogs as described in 
iMagnelli et al.l (|2009[) . 

2.3. Herschel-detected quasars 
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Figure 3. Distribution of the separation between the SDSS po¬ 
sitions and the matched SPIRE 250 pm positions when expanding 
the matching radius to b". The left panel shows the separation 
as a function of the SNR of the 250 pm detection. The solid and 
the dashed lines are the theoretical ~ 1 cr and ~ 2 cr positional un¬ 
certainties {(Jpos)- The circles in yellow color are the sources with 
SNR 250 > 5, while the light blue ones represent the rest. Among 
the yellow circles, the open and the solid ones are those deemed 
by our visual inspection to be “blended” (i.e., affected by source 
blending) and “clean” (free of source blending), respectively. The 
horizontal thick black line going through both panels indicates the 
adopted matching radius of b". 

Our IR quasar sample was derived by matching the 
positions of the SDSS quasars to those of the Herschel 
sources in the band-merged, 250 pm-based catalogs as 
described above. 


2.3.1. Matching radius 

The matching was performed using TOPCAT/STILT^. 
We used a matching radius of 3 ", which is justified below. 

While the Herschel instruments have large beam sizes 
0, the source centroids can still be determined to high 
accuracy. The positional uncertainty of a given source 
depends on its s ignal-to-noise ratio (SNR) (see, e.g., 
Ilvison et al.ll2QQ7[ ). which follows 


Aof = A6 = 


0.6(9 
SNR ’ 


( 1 ) 


where Aa and AS are the nominal 1 a uncertainties of 
RA and Dec, respectively, and 0 is the beam size. The 
SPIRE 250 pm beam size is ^ = 18", which means that 
the 1 a positional uncertainty of a given 250 pm source is 


apos = ^/{Aaf~+jA^ 


15.27" 
SNR ■ 


( 2 ) 


A matching radius of 3 " thus corresponds to ^ 1 cr un¬ 
certainty for o^ects with SNR > 5, or ^ 0.59 cr for those 
with SNR > 20. To further validate our choice, we per¬ 
formed a test using an enlarged matching radius of 5", 
and Figure [ 3 ] demonstrates the results. The separations 
between the 250 pm and the SDSS positions versus the 


^ TOPCAT http://www.starlink.ac.uk/topcat 
STILTS http: //www. starlink. ac. uk/stilts 

5 The FWHM beam sizes are 18", 25" and 36" for the SPIRE 
250, 350 and 500 pm, respectively, and 6-7" and 11-14" for the 

PACS 100 and 160 pm, respective ly. _ 

® We note that ISmith et al.l (|201in use a somewhat differ¬ 
ent relation between the astrometric uncertainty and the SNR, 
(Tpos = 0.655^/SNR, which means that 1 cr uncertainty would be 
2. 4 for SNR = 5 in 250 pm. Our results shown in Figure [3] indicate 
that such a matching radius would be slightly too stringent, and 
therefore we adhered to our choice. 


250 pm SNR are indeed consistent with the expectation 
from Equation ([2]), and the vast majority of the matches 
within 3 " fall within the Icpos curve. Furthermore, the 
separation shows a double-peak feature roughly divided 
at 3 ", indicating that the matches beyond this point are 
likely affected by other factors, such as the source blend- 
i ng problem (s e e belo w). 

IWang et al ] (IMl) have addressed the source posi¬ 
tional accuracy in the HerMES DR2 catalogs through 
end-to-end simulations. For the SF catalogs in 250 pm, 
they find that the real matches between the input and 
the output have the positional offsets peak at around 
5". We note that this is derived using the real matches 
at all SNR levels. If a SNR > 5 threshold is applied, this 
peak is likely to shift to a smaller value. In any case, 
the matching radius of 3" that we adopted is a rather 
conservative choice. 

In a number of fields the PACS data are also available. 
The H-ATLAS band-merged catalog in the SDP field al¬ 
ready includes the PACS photometry (see ^2.2.2^ and 
thus no further action was taken. The HerMES COS¬ 
MOS, GOODS-N, EGS, and Lockman Hole fields have 
PACS data from the PEP program (see 92.2.4p . and we 
matched the coordinates in the PEP 100 and 160 pm 
“blind” catalogs to both the SPIRE “xID250” catalogs 
and the SDSS quasars, again using a matching radius of 
3". 


2.3.2. Source blending 

Due to the large beam sizes of the instruments, Her¬ 
schel images still suffer from a severe source blending 
problem. To evaluate how this problem could impact 
our sample, we visually inspected the Herschel and the 
SDSS images of all the matches in the test case shown in 
Figure O i.e., using a larger matching radius of 5". We 
searched for the signatures of possible source blending, 
such as the presence of close companions in the SDSS im¬ 
ages, and the offset of the source centers among the Her¬ 
schel bands, etc. For this test, we only used the sources 
that have SNR > 5 in 250 pm. The result is also shown 
in Figured As it turns out, enlarging the matching ra¬ 
dius from 3 " to 5 " will include 47 more objects, however 
only 16 of these 47 objects (34%) are “clean” cases. For 
comparison, 83% of the sources within the matching ra¬ 
dius of 3" are “clean”. This lends further support to 
our choice of the matching radius. The drawback of us¬ 
ing this conservative value is that it could exclude some 
objects that have genuine Herschel detections but are 
blended with close neighbor(s). While it is possible to 
in clude such source s after using the methods described 
in lYan et al.l (j2Q14f ) to de-blend, we defer this improve¬ 
ment to our future work. 

In summary, our conservative choice of the matching 
radius has resulted in a sample that is free of significant 
blending problem. This has also simplified the photom¬ 
etry in the SPIRE 350 and 500 pm bands. While the 
source catalogs that we adopted from the survey teams 
were all derived using PSF-fitting based on the 250 pm 
detections, the photometry in 350 and 500 pm would still 
be prone to large errors introduced by blending due to 
their much larger beam sizes (25 " and 36 ", respectively). 
Our relatively clean sample allows us to assume that the 
fluxes measured in these two redder bands are solely con¬ 
tributed by the source detected in 250 pm. 
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T„bb (K) 


Figure 4. Relation between the temperature of the modified 
blackbody spectrum and the temperature inferred from the 

Wein’s displacement law for fho three cases when the emis- 

sivity /3 is 1.0, 1.5 and 2.0, respectively. The dashed line represents 
the equality if these two quantities were the same. 

2.3.3. IR quasar sample summary 

Our Herschel-detected SDSS quasar sample, derived 
using a matching radius of 3" as described above (Ta¬ 
ble [1]), contains 3543 objects in total. Most of our 
studies are based on a subsample of those, called the 
“SNR5” sample, which only includes 134 quasars that 
have SNR > 5 in 250 pm and thus is the sample of the 
most robust Herschel detections. Based on this SNR5 
sample, we further imposed a cut in the 250 pm flux den¬ 
sity, 5'250 ^ 56.6 mJy, and formed a “bright” subsample 
of 82 objects for various discussions below. This flux 
density threshold was adopted based on the lowest flux 
density of the SNR5 objects in the shallowest HerMES 
field L6-XMM-LSS-SWIRE. 

3. DUST EMISSION MODELING 

We inferred the dust emission properties of the IR 
quasars by fitting their EIR SEDs. The thermal emission 
over the full IR regime can be viewed as the collective 
result of all heated dust components of various temper¬ 
atures, and the EIR part is dominated by the coldest 
component. In this paper, we focus on this EIR part in 
the IR quasars and hence our conclusions are pertain¬ 
ing to their cold-dust components. We used two dis¬ 
tinct types of models to fit the EIR SEDs, one being 
a single-temperature, modified blackbody (MBB) spec¬ 
trum and the other being three different sets of star- 
burst templates. Eitting an MBB spectrum to the EIR 
SED is always valid, regardless of the exact dust heating 
sources (i.e., due to photons from either star formation 
or AGN activity). However, it has the drawback that 
the SED must be properly sampled in order to obtain 
well constrained results. Eitting starburst templates, on 
the other hand, is only appropriate if the EIR emission 
is dominated by heating from star formation, and the 
motivation of using these models was to test if the EIR 
emissions are consistent with being caused by star forma¬ 
tion. These two types of SED fitting approach provide 
independent check to each other, and we will show later 
that they also lead to insights into the heating sources. 

3.1. SED fitting using MBB model 


^ We note that there are 18 more matched objects from HerMES 
DR2 catalogs discarded due to SNR250 < 1 


We used the cmcirsed cod43 of [Case 3 (j2Q12f ) to per¬ 
form the MBB fitting, which allowed us to derive the 
IR luminosity, the dust temperature, and the dust mass. 
This procedure was only carried out for the objects that 
have photometry in all the three SPIRE bands (187 ob¬ 
jects in total, among which 103 are in the SNR5 subsam¬ 
ple), because the MBB fit would become unconstrained 
with less bands. 

The EIR emission due to MBB can be written as 


5 (A) 


(l_e-(¥)")(£)3 

Qhc/(\kT^t,t,) _ I 


(3) 


where Tmbb is the characteristic temperature of the MBB, 
Nmbb is the scaling factor that is related to the intrinsic 
luminosity, /3 is the emissivity, and Aq is the reference 
wavelength where the opacity is unity. As most of our 
quasars only have three SPIRE bands available, we had 
to limit the degrees of freedom. We adopted the default 
emissivity of /3 = 1 .5, which is the value typically as¬ 
sumed for cold dust ([Case illoil. By default, cmcirsed 
sets An = 200 inn . We adopted Aq = 100 pm, follow¬ 
ing iDraind (|2006[ ) . While the exact choice of Aq only 
marginally affects the estimates of the total IR luminos¬ 
ity and the dust mass, it will significantly impact the 
estimate of the dust temperature. We will further dis¬ 
cuss this effect in Appendix |Al 
We note that the above form is for general opacity. In 

the optical thin case, at A Aq, the term (1—e“^^^^) re¬ 
duces to (^)^, which is often adopted in the submm/mm 
regime. Throughout this work, we used the general opac¬ 
ity form as in Equation ©• 

The cmcirsed code has the capability of superposing 
a power-law component (PL) to the MBB spectrum to 
accommodate the possible warm dust component whose 
effect could be present in the mid-IR regime (typically at 
< 50 pm in restframe). Thirteen quasars (seven of them 
are in the SNR5 sample) have PACS data in addition to 
the three band SPIRE data, and we utilized this capabil¬ 
ity when fitting these objects. In this case, the MBB-fPL 
model then reads 


^(A) = m, 


Qhc/ (XkT^bb) _ 2 


NpiX 


a V Ar ^ 


(4) 


where a is the PL slope, Npi is the normalization of the 
PL part, and A ^ is the turnover wavelength (for detail, 
see ICasevl[2QT2h . Again, to limit the degrees of freedom, 
a had t o be fixed, a nd we adopted the default value of 
a = 2.0 (ICa,sevll2ni2n . 

Recalling that the EIR emission is dominated by the 
cold-dust component, we can obtain the total IR lumi¬ 
nosity of this component as 

^1000/im 

A)dA (5) 

Jsjum 

by integrating the best-fit model from 8 to 

1000 pm. Eor the 13 objects that have PACS data, we 
also calculated their total IR luminosities as 

^1000/im 

L/fl = = / £'"'''’+P'(A)dA. (6) 

Js/im 

® http://herschel.uci.edu/cmcasey/sedfitting.html 
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Redshift 


Figure 5. Derived IR luminosities of the IR quasars in our sample, where “A”’ denotes one of the four models in use, namely, 

MBB(+PL), SK07, CEOl, and DH02. The errors between the MBB results and those from the starburst templates are not directly 
comparable because of the different approaches adopted in evaluating the errors. The colored symbols represent the SNR5 sample, 
while the grey squares represent the rest. Among the SNR5 objects, the red diamonds indicate those that are in the bright subsample 
(‘S'250 ^ 56.6 mJy), while the yellow circles indicate those that are not. The dark green solid circles in the MBB panel are the objects with 
BAGS data available and hence a power-law component was added to the MBB spectrum in the fitting. To illustrate the impact of the 
survey limit, the limits of corresponding to the fiducial fiux density limit of S'250 = 56.6 mJy are shown as the dashed lines in the three 
panels for the starburst template fits (SK07, CEOl and DH02). In each of these cases, the limit is derived from the entire library by using 
the template with the lowest possible at a given redshift. In the panel for the MBB fit, the limits are given using three different 
of 15, 25, and 35 K, shown as the dashed, the dot-dashed and the dotted lines, respectively. 
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“ 80 
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Figure 6. Distribution of of fho best-fit models for all the IR 
quasars in our sample, using the four different methods as detailed 
in the text. 

We emphasize that 7//^^ (here represented by the MBB 
fit) is the luminosity of the cold-dust component over the 
entire IR range ( 8-1000 pm). While the bulk of its emis¬ 
sion is in the FIR regime, this cold-dust component emits 
beyond FIR and thus it is necessary to integrate over 
8-1000 pm to capture its total IR luminosity. Note that 
this is not the total IR luminosity {Ljr) of the galaxy 
that includes the contributions oJF all dust components 

over 8-1000 pm. In other words, we have < Ljr. 

The best-fit temperature from cmcirsed was taken as 
the temperature of the cold dust, i.e., T^ust = Tmbb- 


We also calculated the “peak” temperature inferred from 
Wien’s displacement law, which is given by 

Tpeak — ^/^peak 5 (7) 

where the coefficient b = 2.898 x 10 ^ pmK. As compared 
to Tjjibb: Tpeak is less sensitive to the specifics of the dust 
emission models in use and therefore could be a better 
proxy to the dust temperature when comparing results 
derived based on different types of templates. For this 
reason, while we mostly use Tdust = T^bb in this paper, 
we also use Tpeak in some occasions. The relation be¬ 
tween Tpeak and Tmbb is shown in Figure IH In case of 
P = 1.5, the relation can be roughly described by the 
following broken linear function: 

_ r0.7653T^66 + 0.9529 , Tmbb < 35 K 
\0.5485T„6b + 8.5153, T„bb>35K. 

3.2. SED fitting using starburst templates 

We also fitted the FIR SEDs using three different li¬ 
braries of starb urst galaxies separately, namely , the theo¬ 
retical model of iSiebenmorgen fc Kriigell (j2QQ7l . hereafter 
SKQ7; 7220 templates ), and the empirical templates of 
iCharv fc ElbazI (l 2 QQlL h ereafter CEOl; 105 templates) 
and iDale fc Heloul (j 2002 l . hereafter DH 02 ; 64 templates). 
Eor a given quasar, the rest frame templates were red- 
shifted according to the quasar’s redshift and convolved 
with the Hersehel passband response curves, and then 
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Figure 7. Comparisons of IR luminosities derived using the MBB model and those based on the three sets of starburst templates. The 
symbols are the same as in Figure^ The upper panels show the comparisons to (he., computed over the full IR range of 8-1000pm), 
while the botto m pa nels show the comparison to (computed over the FIR regime of 60-1000 pm), where “SB” is one of SK07, CEOl 

and DH02. See ^4.II for details. 


were compared to the observed IR quasars SEDs. The 
best-fit template was chosen by minimizing the y-square: 



We note again that the starburst model fitting was also 
done for the objects that have photometry in only two 
or even one SPIRE band. During the fitting, we did 
not re-scale the templates. By construction, these tem¬ 
plates all have associated Ljr values (over 8-1000pm), 
and the value of the best-fit template was adopted as the 
total IR luminosity of the object in question. Eor clarity, 
we will refer to it as where “SB” can be replaced 

by “SK07”, “CEOl” or “DH02” when appropriate. The 
error of Lj^ was estimated by taking the difference be¬ 
tween the Ljr value of the best-fit template and that of 
the one of the second smallest As compared to the 
formal likelihood method, this simple approach has the 
advantage that it works consistently when the param¬ 
eter space is discrete, and that it includes the possible 
systematic errors intrinsic to the template set. 

We should emphasize that Lj^ thus derived is the total 
IR luminosity of the galaxy. If using the starburst mod¬ 
els is appropriate, we should have Lj^ > ^ where 

. While cannot be separated from 
in any of these starburst models, we will show that Lj^ 
can help interpret the origin of T. 

3.3. Dust mass and gas mass 

The MBB fit by cmcirsed also resulted in estimate 
of dust mass (hereafter Mdust)- As this quantity is a 


str ong dependen t of dust temperature {Md oc LirT~^; 
see iCase^ l2Q12f ). it should be used with caution. Eor 
example, in our analysis in 21 we will only use those 
that have Tmbb/^Tmbb > 3. 

Applying a nominal gas-to-dust ratio, the gas mass 
(Mgas) can also be obtained. We adopted the nom- 
inal Milky Way gas-to-dust-mass ratio of 140 (e.g. 
iDraine et al.ll2QQ^ for this work. 

4. RESULTS AND DISCUSSIONS 

The major physical properties obtained in 21 all 
given in the online table accompanying this paper. A 
summary of the information co ntai ned in this table is 
given in Appendix [B] (see Table IbT|). Some examples of 
the SED fitting are also provided in Appendix [0 Here 
we discuss in detail these results, some potential selection 
effects, and their implications. 

4.1. IR luminosity 

Eor clarity, the various flavors of IR luminosities dis¬ 
cussed in 21 summarized here: 

• Ljr: the general designation of the total IR lumi¬ 
nosity over 8-1000 pm; 

• : the contribution of the cold-dust component 
to the total IR luminosity over 8-1000 pm, which 
is also referred to as the total IR luminosity of the 
cold-dust component; 

• the MBB best-fit to the EIR SED (as repre¬ 
sented by the SPIRE data points) integrated over 

8-1000 pm, and by definition 
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Figure 8. No correlation between and the absolute magnitudes of quasars (upper panels), their black hole masses (middle panels), 

or the X-ray luminosities (lower panels). The results are shown in four redshift bins as labeled, and the error bars are omitted for clarity. 


• the MBB+PL best-fit to the FIR SED 
(as represented by the SPIRE and the PACS data 
points) integrated over 8-1000 pm, which is the 
measurement of Ljr using the MBB-j-PL model; 

• Lj^: the measurement of Ljr using the star- 
burst models (“SB” is one of “SK07”, “CEOl” and 
“DH02”, depending on the model set in use), and 
effectively is the best-fit SB template integrated 
over 8-1000 pm. 

Eor one of our purposes later, we also define the corre¬ 
sponding quantities in the EIR instead of over the entire 
IR range, i.e., by integrating over 60-1000 pm only. We 
designate these quantities with the subscript of “EIR”, 



The derived LJ^r^ and Lj^ (referred to as L^r) val¬ 
ues of our IR quasars are shown in Eigure [SI and the 
distributions of the best-fit are shown in Eigure |6l re¬ 
spectively. We note that deriving L'^r (and hence L^jr) 
was not always possible because the MBB fit requires the 
EIR SED being well sampled by the three SPIRE bands, 
while obtaining Lj^ could al ways b e do ne because of 
the nature of the method (see We also note 

that the SNR5 subsample, as expected, has the smallest 
errors in Lfj^. In addition, the majority of the objects 
outside of the SNR5 sample still have < 10 and thus 
are also deemed as having reliable measurements. 

Regardless of the exact model set in use, the majority 
of our IR quasars have obtained good fits and the de¬ 
rived IR luminosities also agree to the extend that we 


expect. This is more clearly shown in the upper panels 
of Eigure 13 where Lj^ are compared against L'f'R • On 
average, L'f'^ values are lower by 0.13, 0.23 and 0.25 dex 
as compared to Lj^ from the SK07, the CEOl and the 
DH02 models, respectively. As explained earlier, this 
is due to the fact that the MBB model only includes 

the cold-dust component {Lf^R = L^ir)-, whereas the 
starburst models also include all other components of 
higher temperatures and thus give the total IR luminos¬ 
ity (Ljr). This is also demonstrated by the 13 objects 
that have PACS data (dark green points), for which we 
carried out the MBB+PL fit and thus obtained the total 
IR luminosity in the form of As one can see, 

there are no offsets between and Lj^. 

To further demonstrate this point, the lower panel of 
Eigure [7] shows similar comparisons as in the upper panel, 
but between L'^j^ and The agreements are excel¬ 

lent and no systematic offsets are found. This can be 
understood as follows. The emission from the hot-dust 
components should be minimal in the EIR regime, and 
hence integrating the starburst models in the EIR should 
only capture the contribution from the cold-dust compo¬ 
nent, i.e., = L^rjr = L^ir- 

As we have been emphasizing, the MBB fit is inde¬ 
pendent of the heating source, while the SB fit being 
valid hinges upon the heating source being star forma¬ 
tion. Therefore, the agreement between L'^j^ and 
strongly suggests that the SB fits are valid and that the 
EIR emissions in these IR quasars are due to star forma- 
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tion. The corollary then is that L'f'^ (= L^j^) is due 
to star formation, because it is the total luminosity of 
the same cold-dust component that gives rise to the FIR 
emission. In the rest of this paper, the quoted IR lumi¬ 
nosity is unless explicitly stated otherwise (mostly 

in ^4.5114.Gl and Appendix [P]) . and we use and 

interchangeably depending on the context. While L'f'^ 
could underestimate the true Lir (see the top panel of 
Figure O by a factor of 1.35 (as compared to to 

1.70-1.78 (as compared to L^r^^ or we prefer 

to be conservative due to the lack of observational con¬ 
straints in the mid-IR for our entire sample. This, how¬ 
ever, is not necessarily a drawback because the mid-IR 
emission, unlike the FIR one, could be seriously contam¬ 
inated by the AGN contribution (see Appendix IC]) . 

Our IR quasars have values ranging from ^ 10^^-^ 
to 10^^'^ Lq (after discarding two objects whose SEDs are 
barely constrained). Most of them (> 80%) are ULIRGs 
{Lir > 10^^ Lq), and some of them (> 15-23%) are even 
HyLIRGs {Ljr > 10^^ Lq). As Figure [5] indicates, there 
is a trend of versus redshifts. Obviously, the lack of 
IR quasars with low at high redshifts is caused by the 
selection effect due to the survey limit. For illustration. 
Figure [5] shows the selection limit corresponding to 
a 250 pm flux density limit of 5'250 = 56.6 mJy (which is 
what we adopted to select the bright subsample from the 
SNR5 sample). Interestingly, there seems to be a deficit 
of very luminous IR quasars at 2 ; < 1, which reflects a 
genuine IR luminosity evolution that is broadly consis¬ 
tent with the evolution of ULIRGs, i.e., there are more 
ULIRGs at 2 ; > 1 than at lower redshifts. 

Our conclusion that the FIR emission of IR quasars 
are powered by star forming activity in dust-rich en¬ 
vironments has also been su ggested by previous stud¬ 
ies at high redshifts (e.g., IWang et al.l l2Qllb[ ). If 
this is indeed the case, usin g the standard Ljr to 
SFR conversion of iKennicnttl (jl998f ). i.e., SFRjr = 
1.0 X Lir/Lq for a Ghabrier initial mass function 
(IMF)0, the SFR of the HyLIRGs in our sample would 
be - 1.0-6.3 X 10^ Mq yr-^ 0. 

There might be concerns whether such extreme SFRs 
are physical. The SFR could indeed be overestimated 
in two ways. First, one could argue that AGN heating 

is still an important contributor to L^^^ and hence the 
SFR cannot be calculated without subtracting this con¬ 
tribution. While there is no viable model quantitatively 
showing that this could be the case (in fact, all available 
models assume the opposite), we cannot yet assert that 
this is impossible. Our argument of L'pfj^ = Lp^^ pre¬ 
sented earlier is only a necessary condition that L^^p is 
due to the heating from star formation but not a suffi¬ 
cient one, and therefore we cannot rule out such a pos¬ 
sibility based on this argument alone. However, we can 
demonstrate that AGN heating is very unlikely dominant 

in L^jp. If it is dominant, it is expected that L^^p should 


® The conversion would be a factor of ~ 1.7 higher if using a 
Salpeter IMF, which was adopted in IKennicnttl (119981 k 

The most conservative SFR estimates would be using 


(over 60 to 1000 pm) instead of (over 8 to 1000 pm) 

would reduce the SFR values by a factor of ~ 1.5. 


, which 


be positively correlated with the AGN activity, i.e., the 
stronger the AGN is, the larger L^j^p should be. Figure [8] 

shows l\^p versus the quasar absolute i-band magnitude 
(nor malized to 2 ; = 2 , adapted from iShen et al.l (j2Qlll ) 
and iParis et al ] (j2Q14f ) for DR7 and DRIO, respectively, 
and all are based on PSF magnitudes after the Galactic 
extinction correction) in four redshift bins. Apparently, 
no such a correlation can be seen. In the lowest redshift 
bin, the distribution of the objects is completely chaotic. 
In the bins at higher and higher redshifts, we tend to see 
only those objects that are more and more IR luminous, 
which is simply due to the selection effect in a flux-limited 
survey. Even among these the most luminous ones, no 

correlation among L^j^p and Mi can be vouched for. Eur- 

thermore, Eigure[8]also shows L^jp versus the back hole 
mass (Mrh ) for the quasars t hat have these estimates 
(taken from IShen et ah ] (|2QllD for the DR7Q quasars). 
Similarly, no correlation exists. Einally, the bottom pan¬ 
els of Eigure [8] show L^^p versus the hard-band X-ray 
luminosity in the restframe 2-10 keV (T^^^ ^®^) for a lim¬ 
ited number of quasars that we can derive this quantity 
based on the data available in the literature 0. As there 
are only 40 such quasars, we split them into two redshift 
bins, 0 < ^ < 1.7 and 2 ; > 1.7, respectively, so that each 
bin receives approximately the same number of objects 
for statistics (19 and 21 objects, respectively). They all 
have > 10^^ erg/s, which is well above the con¬ 

ventional X-ray AGN selection threshold of 10^^ erg/s, 
above which the X-ray luminosity is believed to be pre¬ 
dominantly due to AGN. Therefore, is a strong 

indicator of the AGN activity. Again, no correlation be¬ 
tween L^j^p and can be seen. This is also very 

consis tent with the rec e nt res ults of iSvmeonidis et all 
(12011) and lAzadi et al.l (I 20 TI ) in the similar 
regime. 

Therefore, while we do not have direct evidence to as¬ 
sert that AGN has no contribution to we do have 

evidence (albeit still indirect) against that AGN con¬ 
tribution can be dominant. This further strengths our 
conclusion of L^j^p being due to star formation based 
on the earlier argument of L^/p = Lp^p. As an addi¬ 
tional check of consistency in our conclusion, we have also 
tested the AGN/s tarbnrst decompositio n approach, us¬ 
ing the method of iMnllanev et all (j2Qllf ) on the objects 
that have data in the PAGS and/or the Spitzer MIPS 
bands. The results are detailed in Appendix O Note 
that all the AGN/starbnrst decompositio n schemes avail¬ 
able i n the literature to date (including IMnllanev et al.l 
(j2Qllf )) assume that the AGN contribution drops off in 
the EIR, and hence the decomposition is not entirely ap¬ 
propriate in asserting whether AGN contribute signifi¬ 
cantly to Lrir. Nevertheless, our result shows that the 
starburst-contributed IR luminosities as derived in the 

were derived based on the data from the Chandra 
Source Catalog Release 1 (IE vans et al.l[207ol) and the 3XMM-DR5 
catalog (|Roseii et al.l[207^ Briefly, a power-law in the form of 
oc 1/“*^ was fit to the flux densities at different energy bands, 
and the total energy in restframe 2-10 keV was calculated by inte¬ 
grating the best-fit power-law over this energy range. The best-fit 
a has a median of ~ 0.7, which correspond to the photon index 
F - 1.7. 
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Figure 9. Distribution of the derived dust temperatures 
Only objects with are shown (134 in total). The 

colored symbols (86 objects) represent those that are within the 
SNR5 sample, while the grey squares (48 objects) represent the 
rest. Among the SNR5 objects, the red diamonds (55 objects) 
indicate those that are in the bright subsample (S'250 > 56.6 mJy), 
while the yellow circles (31 objects) indicate those that are not. 
The right panel shows the histograms of for the displayed 

objects (all: grey; SNR5: yellow; “bright” SNR5: red). 

decomposition scheme (designated as are also 

consistent with our and Lj^ as derived using the 

SPIRE data alone, and therefore we do not find any ev¬ 
idence against our conclusion. 

The other possibility is that the most luminous ob¬ 
jects are actually gravitationally lensed, which means 
that their intrinsic luminosities must be lower and so are 
their SFR estimates. Currently, we do not have further 
data to address this issue. In ^4 31 however, we will show 
that it is also unlikely that the most luminous objects 
are predominantly the result of lensing. 

4.2. Dust temperature 

Due to the limitation of our data, it was not always 
possible to obtain a good temperature estimate using ei¬ 
ther Tmhh 01 * Tpeak- To constraiu the temperature well, 
at least the three SPIRE bands (250, 350 and 500 pm) 
should sample the peak region of the EIR SED on both 
the “rising” and the “falling” sides (i.e., short and long 
wavelength sides, respectively). However, this is not al¬ 
ways satisfied in our cases. To secure our discussion in 
this section, here we only include the objects that have 
Tmbb/^Tmbb ^ 3 (i.c., equivalent to “SNR” in tempera¬ 
ture estimate > 3). 

The left panel of Figure [9] shows the distribution of 
Tmbb from the MBB fit. The majority of them are in the 
range of 30-50 K, which is generally consistent with the 
picture that the FIR emission in these objects is domi¬ 
nated by star formation activity heating the dust. There 
also seems to be an “evolutionary trend” in that Tmbb 
increases with redshifts, however this is due to the bias 
in how a secure estimate of temperature could be de¬ 
rived in our case: the peak of the MBB FIR emission at 
the same temperature would shift to longer wavelengths 
with increasing redshifts, and therefore to keep the peak 
being well sampled by the three SPIRE bands on both 
the rising and the falling side of the SED, a higher tem¬ 
perature would be needed to shift the peak to shorter 
wavelengths. 

We carried out an extensive simulation using the MBB 
models to further investigate this temperature bias, and 
the results are summarized in Figure [TOl To connect to 


the argument above, Tpeak is used for this demonstra¬ 
tion. For a given our goal is to study the con¬ 

straint on Tpeak when the detection limit is imposed. We 

generated a large set of MBB models, with rang¬ 

ing from 10 ^^ to 10 ^^ Lq (in a step-size of 0.2 dex). At 

each we varied Tmbb from 10 to 100 K (step-size 

IK). These models were redshifted to 2 : = 0 to 6 (step- 
size O.I), and the flux densities 5'250 were calculated from 
these simulated spectra. We then imposed a desired de¬ 
tection limit, selected the simulated objects that have 
-5250 above this threshold, and, after calculating their 
Tpeak ^ pat them on the Tpeak — z plane as shown in Fig¬ 
ure nni For illustration purpose, here we use a fiducial 
threshold of ^250 = 56.6 mJy, and we only show three 

cases in namely, = 10^^, 10^^ and 10^^ Lq, 

respectively. The area filled in pink shows the region 

occupied by the = 10 ^^ Lq objects thus selected. 

For the same detection threshold, simulated objects of a 
smaller luminosity will occupy a smaller region. For ex¬ 
ample, the green region within the pink region is where 

the objects with = 10 ^^ Lq reside, and the blue 

region within the green region is where the objects with 

~ ^0 i* 6 side. We refer to such a color-coded 

region as the “occupation region” of a given (Tj^\ ^ 250 ) 
pair. These regions have boundaries at both high and low 
Tpeak ^ 01 * in other words, for the adopted 525 o, an object 

of the given can be detected at a specified red- 

shift only when its dust temperature is within the shown 
boundaries. Lowering 5250 , or equivalently, increasing 

would expand the boundaries of an occupation re¬ 
gion in the way that the blue region would “grow” to the 
green region and to the pink region. 

Note that the occupation regions have distinct 
“bumps” whose tips are aligned on a straight line, which 
is shown in red in Figure [TOl This “ridge” line indicates 
the direction towards which an occupation region would 
“grow” fastest in area when decreasing 5^50 and/or in¬ 
creasing Changing and/or 5^50 will only shift 

the bumps along the ridge line, and the ridge line itself 
does not change as long as the bands involved in mea¬ 
suring the dust temperature stay the same. 

To improve the above simulation further, we added 
photometric errors (according to the HerMES survey lim¬ 
its) to the synthesized SEDs derived from the simulated 
spectra, and fit them using the MBB models. As it turns 
out, the simulated objects that have reliable tempera¬ 
ture estimates {Tmbb/^Tmbb > 3) through the fit are all 
distributed on the narrow ridge line that connects the 
“bumps”. This is shown in the right panel of Figure [9] 
for the case of SNR = 5 in 250 pm. This again shows 
that the dust temperature estimates based on the SPIRE 
data prefer some certain temperature at a given redshift, 
although the SPIRE deteetions can span a wide temper¬ 
ature range. Using a formalism similar to Equation ([3), 
this “preferred” temperature as a function of redshift can 
be approximately described by the following equation: 


^pref 2.898 X lO^pmK 
~ 303.7pm 


( 10 ) 


where 303.7 pm is the “preferred” peak wavelength due 
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Figure 10. Demonstration of the redshift bias in the dust temperature derived using the SPIRE bands. Tpea/c is used for this purpose, 
and for simplicity a fiducial sensitivity threshold of = 56.6 mjy is adopted in this figure. Using the MBB fit of the cmcirsed tool, 
Tpeak can be directly derived, and the results for the SNR5 sample are shown, with those from the bright subsample (S'250 > 56.6 mJy) 
being coded in red diamonds and the rest being coded in yellow dots. The adopted 5250 imposes a selection boundary in the (Tpea/c? 
plane for a given limit, and thus forms a “occupation region” for the objects whose are less than this limit. The shaded color 


regions denote such regions for the limits of 10^^ (light blue), 10^^ (light green) and 10^^ Lq (pink), respectively. Such occupation 

regions have “bumps” whose tips are aligned on a straight ridge line (shown in red), and this ridge line indicates the temperatures preferred 
by the SPIRE bands at different redshifts. The right panel shows the results from the simulation as detailed in the text. The blue curves 
are the number density contours of the simulated SNR = 5 (in 250 pm) objects that have reliable Tpea/c estimates through the MBB fitting, 
and the most dense contours are all distributed close to the ridge line, which further confirms the redshift bias in the dust temperature 
derived based on the SPIRE bands. 


to the usage of the three SPIRE bands. In other words, 
this means that the closer the peak of the FIR emission 
is to 303.7pm/(1 + z), the more reliable the temperature 
estimation will be. 

4.3. Relation between dust temperature and IR 
luminosity 

Here we examine the T^ust-L^iR relation of IR quasars. 
For simplicity, we use in this discussion, and only 
include the 134 objects that have good estimates of both 

< 10) Tmbb (“SNR”> 3). This is shown 
in Figure im where the crosses are the individual ob¬ 
jects and the red circles represent the average at a given 

(step-size of 0.2 dex in W e also plot 

the mean result from iSvmeonidis et al.1 (j2013L light green 
circles), who have analyzed a sample of IR luminous 
{Lir > 10^^ Lq) galaxies at 0.1 < 2 ; < 2 using the deep 
PACS and SPIRE data in the COSMOS, the GOO DS-N 
and the GOODS-S fields. iSvmeonidis et al.l (j2013f ) find 
that their Ljr-T relatiorQ has only a modest increasing 
trend towards high luminosities, and their interpretation 
is that the increase of Ljr is caused by an increase in the 
dust mass and/or the IR emitting radius rather than by 
an increase in the intensity of the dust heating radiation 
field. 

However, the picture is rather different for our sam¬ 
ple because our IR quasars span a much wider range in 

luminosity. Our r e lation agrees reasonably 

with that of ISvmeonidis et ^ (j2013f ) in the overlapped, 

low luminosity range {L^j^r ^ 10^^ L©), and then dra- 

We note that ISvmeonidis et ahl (1201311 adopt Aq = 100 pm 
and /3 = 1.5 as we do here. 


matically rises to higher luminosities. The existence of 
such a relation is against the possibility that our sample 
could be significantly affected by gravitational leasing, 
because the magnificat ion c annot be correlated with the 
dust temperature (see ^4.ip . 

We argue that this relation cannot be at¬ 

tributed to the selection effect of our sample. To demon¬ 
strate this point, we simulated a large number of objects 

of different T^bb and over the redshift range of our 
sample, and recovered them using various selection crite¬ 
ria in Tmbb and 5'250- Figure [12] shows the results in four 
redshift bins, for two different Tmbb “SNR” thresholds of 
three and five, and two different 5'250 thresholds of 1 and 
50 mJy, respectively. The key points can be summarized 
as follows. First, adopting a higher Tmbb “SNR” (e.g., 
Tmbb > ^T^bb instead of Tmbb > ST^^b) would be against 
objects with high Tmbb- Second, adopting a higher 5'250 
threshold (e.g., 5'250 > 50mJy instead of 5'250 > ImJy) 

would be against objects with low L^ir- To reiterate, 
our current work adopts Tmbb > ‘^'^mbb- While our sam¬ 
ple does not have an uniform 5'250 threshold due to the 
varying survey limits in different fields, our objects all 
have 5'250 > 50mJy. From Figure O one can see that 
the objects that have highest probability of being selected 

by our criteria would be those with higher and lower 
Tmbb than our data points, however such objects are not 
presented in our sample, i.e., there is a genuine lack of 
such objects. 

Before discussing the lack of objects with (high 
low Tmbb)i let us first understand the increasing trend of 

Tmbb with increasing T. Recall that for a perfect black 
body, Stefan-Boltzmann law states that L = AirR'^aT^. 
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Figure 11. Strong relation as inferred from our IR quasars that have reliable and measurements. The light grey 

crosses represent the individual quasars (134 objects in total), while the big red circles show their mean values in bins of 0.2 dex. 

The thick error bars on these red circles are the standard deviations around the mean values, while the t hin error bars are the summed (in 
quadrature) errors of all objects in each bin. For comparison, the results from (|Svmeonidis et aI1l2013l j are shown as the green symbols. 
The dashed lines represent the MBB equivalent of Stefan-Boltzmann law as derived based on Equation m, using a family of R^f f values 
as labeled. 


Motivated by this, we integrated Equation (j3j) and found 
that the equivalent for a general opacity MBB should fol¬ 
low , where Reff of a given galaxy 

should be interpreted as the effective radius of the equiv¬ 
alent FIR emitting region if we combine together all its 
dust-enshrouded star-forming regions. As shown in Fig¬ 
ure El our data points can be explained by this relation 

with a family of Reff- At < 10 ^^ Lq, the increasing 

in is mostly dominated by the increasing in Reff, 

which range from ^ 0.1 to 0.5 kpc, and the increasing in 
Tjjibb only plays a modest role. This is consistent with 
the suggestion of iSvmeonidis et ahl (j2013[ ) as summarized 
earlier. Naturally, Reff cannot be increased indefinitely 
because the sizes of galaxies are finite. Our data suggest 
that Reff reaches its maximum of ^ 0.5-0.7 kpc at the 

ULIRG luminosity. At > 10^^ Lq, the increasing in 

taken over by the increasing of Tmbb, which can 
be due to more intense radiation field caused by more 

intense starburst activity. The lack of (high low 

Tmbb) objects thus is the result of the limit in Reff- 

This also suggests that the Tj^bb-L^iR relation as seen 
in Figure [TT] for IR quasars is an envelope of the general 

distribution of IR-luminous objects on the {Tmbb, L^ir) 

plane. In other words, for a given L^jr, the dust tem¬ 
perature reached in IR quasars is the lowest among all 
possibilit ies in IR-luminou s obje cts. In fact, the data 
points of iSvmeonidis et ^ (j2Q13[ ) indeed are above ours 
in Figure dlj which is perfectly consistent with our inter¬ 
pretation. However, we do not have an explanation on 
why this envelope manifests itself in IR quasars. 


4.4. Dust mass and gas mass 

The MBB fits also resulted in estimates of dust mass 
(hereafter Mdust), whose distribution is shown in Fig¬ 
ure [13] with respect to and Tmbb- As the calculation 

of Mdust is stron gly affected by Tmbb {Mdust oc 
see iCasevI |2Q12[) , again only those with T^bb/^mbb — ^ 
are included in the plot. The distribution of Mdust peaks 
at ^ Mq, which seems to be higher than the 

dust cont ents of ULIRGs at these redshifts in general 
(see, e.g. JYan et al.ll2Q14l) . However, considering that a 

good fraction of our objects have L^j^r > 10^^ Lq, such 
high dust masses probably are not surprising. Interest¬ 
ingly, for the objects in our sample. Figure [13] also sug¬ 
gests that Mdust is almost a flat distribution of Tmbb- 
This can be explained by the observed T^bb-L^iR rela¬ 
tion of our IR quasars, which can be approximated by 

^ where a ^ 4-5. This means that for this 

group of objects we should observe a flat distribution of 
Mdust with respect to Tmbb, which is exactly what Fig¬ 
ure shows. Therefore, our results are self-consistent. 

Adopting a nominal gas-to-dust ratio of 140, we ob¬ 
tain the gas masses of these objects, which peak at 
Mgas ^ (0.8-4.4) X IO^^Mq. This indicates that the IR 
quasar host galaxies are very rich in gas. If they could 
turn all their gas into stars in their current ULIRG phase, 
they would grow substantially in stellar masses. In fact, 
the added stars alone would amount to the stellar masses 
of typical giant elliptical galaxies in the local universe, 
which are to the order of 10^^ Mq. Figure [14] shows the 
time scale, tdep, that their host galaxies would deplete 
the gas reservoir if they would keep forming stars at the 
rates as seen in their current ULIRG phase. Most of these 
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Figure 12. Simulation results showing that the observed relation is not due to selection effect. As labeled in each panel, the 

results are given in four redshift bins (in column) for the combinations of two thresholds and two 5250 thresholds (in row). The red 

symbols are the same relation as in Figure [TTl while the blue blocks in the background represent the densities of the simulated 

objects recovered by the labeled selection thresholds (the coding of the color depth is shown to the right). Under our current sample 
selection (the first row), the highest densities of the simulated objects all occur below the relation, i.e., in the regions of higher 


and lower than those defined by the 


(cd) 

IR 


In other words, while the objects of higher and lower 


relation. 

would have the highest probability of being included by our selection, they do not present in our sample, i.e., there is a genuine lack of such 
objects in the universe. On the other hand, our selection (mostly the 5250 threshold due to the survey limit) is against the objects of higher 

Tmbb- Therefore, it is most likely that the T^bb-L^j^R relation, discovered among IR quasars, is the envelope of the general distribution of 
IR galaxies. 


objects have tdep ^ 560 Myr (Log{tdep) ^ 8.75), which 
are broadly consistent with the duration of ULIRGs and 
therefore would suggest that they could indeed turn all 
their gas reservoir in their current ULIRG phase. How¬ 
ever, there are a few objects (22, among which 5 are 
among those of the most secure Tmbb estimates) that have 
tdep > 1 Gyr (among which one has tdep > 5Gyr). It is 
unclear whether such objects would be able to keep their 
extreme SFRs over such a long period. 

4.5. Fraction of IR-luminous quasars 

A question of general interest is how many optically 
selected quasars are FIR-bright. Table [T] has provided 
a rough answer to this question: there are 6710 SDSS 
quasars in our Herschel fields, and 354 (5.3%) of them 
have 250 pm detections, among which 134 (2.0%) are de¬ 
tected at SNR >50 However, little can be further 
inferred from such statistics, because this is set by their 
detections in the 250 pm images, which are entangled by 
the varying survey limits over the fields, the cosmolog¬ 
ical dimming effects and the /^-corrections. To address 

Our 250 p m detection rate of 5.3% is lower than that of 
IDai et ahl (I2012|j . who obtain a ~ 10% 250 pm detection rate among 
their 326 quasars in the HerMES Lockman Hole field. Such a dif¬ 
ference is due to the difference in the parent quasar samples. Most 
importantly, the sample of Dai et al. is pre-selected using the 
MIPS 24 pm detections, which presumably is biased to favor FIR 
detections. 


this question in a better defined manner, one should dis¬ 
cuss it in terms of luminosity. To take advantage of all 
the 250 pm detections for this analysis, here we use Ljr 
because its derivation does not require all three SPIRE 
bands (as oppose to L'f'R)- This also means that we 

are discussing in terms of Lir rather than as we 

do in the previous sections. To minimize the systematic 
effects introduced by a particular set of templates, we 
adopt the average of the three derived values based on 
the three sets (SK07, GEOl and DH02). Eor the sake of 
simplicity, we still denote this quantity as Ljr. 

We introduce a “critical” luminosity, Lc, such that 
quasars with Ljr > Lc are deemed as ^RR-luminous^\ 
The exact choice of Lc is somewhat arbitrary, and here we 
adopt Lc = Lq. The redshift distribution of these 

IR-luminous quasars and their fraction among the SDSS 
quasars is shown in Eigure [151 The yellow symbols are 
for the case of the SNR5 sample (i.e., the IR-luminous 
quasars among the SNR5 sample), while the blue and 
red symbols are for the cases of the whole IR quasar 
sample and the bright SNR5 sample, respectively, which 
represent the most aggressive and the most conservative 
inclusions of objects in the calculation. Obviously, the 
fraction of IR-luminous quasars is not flat over all red- 
shifts, which is suggestive of evolution with time. To 
be specific, this fraction peaks at around z ^ 1-2, and 
drops sharply from z ^ 2 to higher redshifts. While the 
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Figure 13. Distribution of dust mass {Mdust) of fho IR quasars with respect to (left panel) and T^bb (right panel). The middle panel 

shows the histogram. The gray squares, yellow circles, and red diamonds are for the objects with 5 > > 3, 10 > > 5, 

and > 10, respectively. The gray, yellow, and red histograms in the middle panel are for the objects with > 3, 5, 

and 10, respectively. 
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Figure 14. Distribution of the gas depletion time t^ep for the 
IR quasars. The left panel shows the distribution with respect to 
redshifts and look-back time. A number of objects have tdep larger 
than the look-back time, which could mean that they would still 
have plenty of gas left when evolve to today. The middle panel 
shows the histogram of t^ep- The symbols are the same as in 
Figure [131 

“spiky” features at 2 ; > 3.5 could be attributed to the 
small numbers in both the IR quasar sample and its par¬ 
ent SDSS quasar sample, it would be difficult to explain 
the notable drop from 2 ; ^ 2 to 2 ; ^ 3 using this factor 
alone. 

Two factors could impact the above picture, however. 
First of all, the exact fraction of IR-luminous quasars of 
course depends on the adopted Lc. From Figure O it is 
clear that increasing Lc would lower the current peak and 
shift it slightly to higher redshifts, because there would 
be fewer objects qualified as IR-luminous at lower red¬ 
shifts. However, Lc should not be increased arbitrarily. 
For example, if we were to adopt Lc = 10^^ Lq, only a few 
objects would have survived and this analysis would not 
be meaningful. On the other hand, lowering Lc would 
not change the fractions at 2 ; > 1 because the number 
of IR-luminous quasars would not change in our sample. 
However, this is largely due to the second factor, namely, 
the survey incompleteness. In fact, again from Figure [5] 


one can see that the incompleteness at Lq could 

become significant at 2 ; > 2.5. We investigated the ef¬ 
fect of this incompleteness by stacking the SDSS quasars 
that did not find matches in the Herschel catalogs (here¬ 
after Herschel-undetected quasars”). The details of this 
stacking analysis is presented in Appendix [D1 
Basically, we need to take into account the objects that 
have Lj^ ^ Lc = Lq and yet are missing from our 

sample due to the survey limits. To focus on the drop 
at 2 ; ~ 2, we only discuss the redshift ranges 1 < 2 ; < 2 
and 2 < 2 ; < 3. While the correction factors cannot be 
determined precisely, we obtained the maximum factor 
possible at 1 < 2 ; < 2, fcor < 5.4, and the minimum 
factor necessary at 2 < 2 ; < 3, fcor >3.1. This could 
potentially increase the fraction of IR-luminous quasars 
from ^ 4.5% to < 24% at 1 < ^ < 2, and will definitely 
increase it from ^ 2.5% to < 7.7% at 2 < z < 3. There¬ 
fore, it is possible that the drop at 2 ; ^ 2 still persists after 
the incompleteness correction, although our current data 
do not allow us determine whether the drop is as steep 
as what inferred from using only the Herschel-detected 
objects. 

4.6. Contribution to luminosity density 
Here we investigate the contribution of IR quasars to 
the IR luminosity density This is done by 

adding Lj^ of the IR quasars in a given redshift bin and 
then dividing by the total survey volume in this bin. We 
do not intend to correct for the incompleteness imposed 
by the Herschel survey limits, and thus what we can ob¬ 
tain would only be a strict lower limit of the contribution 
from the optical quasar population. 

While our sample is the largest one possible at this 
stage, its size is still rather limited when being divided 
into subsamples by redshifts, and thus the exact step- 
sizes in the redshift domain will slightly affect the de¬ 
tailed results. For this reason, we adopted two types 
of binning in redshifts, one being a uniform division at 
a step-size of Az = 0.5 and the other being equal vol¬ 
ume (1.46 X lO^Mpc^) in the successive bins. The re- 
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Figure 15. Evolution of “IR bright” {LjS > Lc) quasars with respect to redshifts (lower panel). In this work, we adopt Lc = 10 12-5 Lq. 
For comparison, the redshift distribution of the parent SDSS quasar sample is reproduced in the upper panel. The colored histograms in the 
lower panel show the distributions of IR bright quasars among the entire sample (blue), the SNR5 sample (yellow), and its bright subsample 
(red). The error bars denote the uncertainties based on the Poissonian statistics. In all cases, there is a peak at around 1-5 < 2; < 2.5, 
which does not follow the distribution of the parent SDSS quasar sample. To demonstrate it more clearly, the curves (with the same color 
coding as the histograms) show the fraction of the IR bright quasars among the parent quasars. The shaded regions denotes the Poissonian 
uncertainties. The fraction has a prominent drop at around z ^ 2. The two arrows indicate the impa ct of the incompleten ess correction as 
detailed in the text and Appendix [D] For reference, such a fraction at z 6 is also calculated using [Leipski et HI (120141) and is shown as 
the yellow asterisk. 
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Figure 16. Contribution of IR quasars to the IR luminosity density as a function of redshift, for the cases of equal redshift binning 
(left) and equal volume binning (right). The IR luminosities used here are based on the starburst models (i.e., see ESI for details). 

The error bars take into account both the errors in the derived values and the Poissonian uncertainties in the number counts. The 

yellow asterisk denotes the result derived using the high-2; IR quasar sample of ILeipski et al.l (120141) . The horizontal dashed lines are the 
contribution of the Herschel-undetected quasars based on the stacking analysis as detailed in Appendix m 


suits are shown in Figure [161 While the detailed features 
are somewhat different in these two schemes, the overall 
characteristics are the same. 

First of all, IR quasars among optical quasar popula¬ 
tion only contribute a very small fraction to the total 


IR luminosity de nsity. One r e cent e xample to compare 
to is the study of IViero et al.l (|2Q13[) . where the authors 
use a K-band selected galaxy sample to perform a simi¬ 
lar analysis in the FIR and find that the IR luminosity 
density produced by galaxies peaks around z = 1-2 at 
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^ 10 ^-^ Lq Mpc“^. As Figure [T6l shows, dust emission 
due to star forming activity of IR quasars only contribute 
^ 0.06% of this amount. 

Second, declines at 2 ; > 2. This echoes the 

decline of IR luminous quasar fraction at the same red- 
shift as discussed above (see Figure [T5]) . However, this 
picture is also affected by the survey limits. To take 
the Herschel-undetected quasars into account, we used 
the stacking results as described in Appendix [Dl While 
we were only able to make corrections over four redshift 
bins with much larger step-size {Az = 1), it is clear that 
after the correction peaks at 2 ; ^ 2.5 instead of 
2 ; ~ 1.5 when uncorrected. 

It is well known that optical qua sar number density 
peaks at around 2 ; ~ 2-3 (see e.g., lOsme rllHol for re¬ 
view). Using the SDSS DR3 quasars. iRichards et al.l 
(|2QQ6[ ) find the peak at 2 ; ^ 2.5 when integrating the 
Tband lum i nosity function to Mi[z = 2] = —27.6mag. 
iJiang et al.l (|2QQ6[) use a faint quasar sample in an area of 
3.9 deg^ within the SDSS Stripe 82 and find that the peak 
would sh ift to 2 ; ^ 2 if in tegrating to Mg = —22.5 mag 
(see also iRoss et al.|[2Q13F The UV luminosity density 
(and hence the global SFR density) of normal galax¬ 
ies rises sharply from 2 ; = 0 to higher redshifts, also 
pea ks at 2 ; ^ 2, but keeps fiat out to 2 ; ^ 4-5 (see 
e.g.. lHoDkins fc Beacom|[2QQ^ . It is thus intriguing that 
peaks at around the same redshift, which sug¬ 
gests that this could be the result of the co-evolution of 
extreme star formation and quasars. 

4.7. Two quasars with increasing FIR SEDs 

In the course of the FIR SED analysis, we found two 
unique quasars that have monolithically increasing SEDs 
from 250 to 500 pnO. They are J090910.08+012135.6 
(“HATLAS-SDP-001” for short) at 2 ; = 1.02 and 

J012528.83-000555.9 (“HerS-080”) at z = 1.08, re¬ 
spectively. By searching the NASA/IPAC Extragalac- 
tic Database (NED) , we found that bot h of them are 
blazars cataloged bv iMassaro et al ] (|2009[ ). with the ob¬ 
ject names of BZQJ0909+0121 and BZQJ0125-0005, re¬ 
spectively. The EIR pr operties of HATLAS-SDP-OOl has 
also been discussed bv iGonzalez-Nuevo et al.l (j2010[ ). 

The monolithically increasing SEDs of these objects 
mimic the SEDs of the so-called “500 pm peakers”, which 
are c andidates of EIR galaxies at 2 ; > 4 (jPope fc CharvI 
IMH). Apparently, such objects can be non-negligible 
contaminators to such EIR high-z candidates, as they 
pass the usual crite rion of S'soo/^'sso > 1-3 (see e.g., 
iRiechers et al.l[20T^ . 

5. CONCLUSION AND SUMMARY 

In this work, we combined the SDSS quasars from 
DR7 and DRIO and searched for their EIR counterparts 
in the 172 deg^ Herschel wide survey fields where the 
high-level Herschel maps have been made public by the 
relevant survey teams. Erom the total of 6170 SDSS 
quasars within the Herschel field coverage, we assembled 
a sample of 354 quasars that are detected in the Herschel 
SPIRE data, 134 of which are highly secure, SNR > 5 
sources in 250 pm. As we used a stringent matching cri- 

These two objects have been excluded from our analysis else¬ 
where in the paper. 


terion, the contamination due to source blending is min¬ 
imal. This IR quasar sample spans a wide redshift range 
of 0.14 < 2 ; < 4.7, and is the largest sample of optical 
quasars that have EIR detections. 

To investigate their properties, we analyzed their EIR 
SEDs using a modified black body model (MBB) as well 
as three different sets of starburst (SB) templates (SK07, 
CEOl and DH02). By focusing on the EIR emission, we 
confine our discussion mainly to the cold-dust component 
of IR quasars, as the EIR emission is predominantly due 
to this component. Our conclusions are summarized be¬ 
low. 

• The results based on the MBB model (independent 
of the heating source) are consistent with those 
based on the SB models (legitimate only when 
the heating source is from star formation), which 
strongly suggests that the IR luminosity of the 

cold-dust component in these IR quasars 
is mainly due to the heating of star formation 
rather than AGN. This is further strengthened 
by the additional supporting evidence that there is 

no positive correlation between and the ab¬ 

solute magnitudes or the black hole masses of the 
IR quasars, both of the latter being indicators of 

the strength of AGN. derived in this work 

as based on the SPIRE photometry, underes¬ 
timates the total IR luminosity because it does not 
include the contribution from other warmer dust 
components. However, as it is very unlikely be¬ 
ing significantly contaminated by the AGN heat¬ 
ing, this quantity is more preferable in inferring 
the SER of the host galaxies. 

• The derived values, adopting the more con¬ 

servative ones based on the MBB fitting results, 
range from to Lq (after discarding two 

objects whose SEDs are barely constrained), with 
^ 80% being ULIRG (> 10^^ Lq) and ^ 15% being 
HyLIRG (> 10^^ Lq). There is a general trend that 

increases at larger redshifts, which is mostly 
due to the selection effect in our flux-limited sam¬ 
ple. However, there is a lack of high objects at 
the lowest redshifts {z < I), which is broadly con¬ 
sistent with the picture that ULIRGs are scarce in 
the low-redshift universe. 

• Due to the wavelength coverage of the SPIRE pass- 
bands, the dust temperature can only be well con¬ 
strained for a fraction of the IR quasars whose 
SPIRE detections sample both the blue and the 
red sides to the peak of their EIR emissions. Eor 
these objects, the derived temperatures show an in¬ 
creasing trend with redshift, which is caused by this 
selection effect. We show that these objects dis¬ 
tribute along a “preferred” region on the dust tem¬ 
perate versus redshift plane, and this temperature 
bias is common to any Herschel sources that are 
limited by the SPIRE data. Nevertheless, we find 
that most of the IR quasars with well constrained 
temperatures have Tmbb ^ 25 to 50 K, which is con¬ 
sistent with the picture that their EIR emission is 
mostly due to young stars heating the cold ISM. 
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Figure Al. Impact of the different choices of Aq to the dust temperature and the IR luminosity estimates, using a representative case 
where the input MBB models using Aq = 100 pm are at 2; = 2 and have = 10^^ L©. The left panel shows the comparison of the output 

Tmhh (blue symbols) and (red symbols) estimates, derived using the MBB fit with Aq = 200 pm, to the input temperature. Similarly, 

the right panel shows the comparison of estimates to the input value. In both panels, the shaded region indicates the temperature 

range within which the peak region of the MBB spectra are well sampled by the SPIRE bands. 


• In spite of the dust temperature bias, the IR 
quasars with well constrained dust temperatures 

allow us to investigate the T-L^j^ relation over 

a wide dynamic range in We find that 

there is a dramatic increase of dust temperature 
at > 10^^ Lq. Through simulations, we show 
that this trend cannot be due to the selection effect 
of our sample. Instead, this trend, which holds for 
IR quasars, seems to be the envelope of the general 

distribution of IR objects on the (T, L^j^) plane. 
At the low luminosity end along this envelope, the 

increasing of is largely due to the increasing in 
the effective radius of the heated region (or equiv¬ 
alently, the enclosed dust mass). The behavior of 
the trend shows that the size of the heated region 
cannot be arbitrarily increased, and any further in¬ 
creasing of must be largely driven by the in¬ 
creased heating (i.e., more intense star formation 
rate per unit volume). 

• The SFR values inferred from range from 6.3 
to 63IOM0yr“^ (for a Chabrier IMF; the values 
would be I.7x higher if using a Salpeter IMF). 
From the dust mass derived via MBB fitting, and 
using a nominal gas-to-dust ratio of 140, we have 
inferred the gas mass for the IR quasars, most of 
which being within (0.8-4.4) X lO^^M©. Given 
the SFR and the gas mass, for most of the objects in 
our sample, the time scale that their host galaxies 
would deplete the gas reservoir is tdep ^ 560 Myr, 
while a few of them could have tdep > I Gyr. 


• The fraction of IR-luminous optical quasars evolves 
with time. For a fiducial threshold of Lc = 
10^^-^ Lq, the fraction of quasars with Ljr > Lc 
peaks at around z ^ 2, an epoch in general agree¬ 
ment with the peak of the global star formation 
rate density evolution. 

• IR quasars, even counting those that are not de¬ 
tected in the current Herschel surveys, only con¬ 
tribute a very small fraction to the total IR lu¬ 
minosity density. This contribution also peaks at 
around 2 ; ^ 2-3. 

The full catalog of our sample, including all the derived 
physical properties, is available as the online data of this 
paper. 
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APPENDIX 


A. IMPACT OF Ao IN MBB MODELS 

As detailed in EH we adopted the MBB model of the cmcirsed code (jCasevll2012tl . which takes the following form: 


5(A) = N^bb 


qHc/ (^XkT^hh) — t 


The term 1 — modifies the black body term in that it includes a wavelength-dependent optical depth r(A), 

which is assumed to follow a power-law of r(A) = (^)^- The parameter Aq is the wavelength where the optical depth 
is unity. By default, cmcirsed sets Aq = 200 pm. In this work, we followed IPraind (j2006f ) and adopted Aq = 100 pm. 

The difference in the choice of Aq impacts the derived T^bb significantly. This is because the “modifying” term 
varies significantly when different Aq is adopted. As an example, let us consider a FIR SED that peaks at restframe 
A = 100 pm. The “modifying” term is 0.63 at this peak wavelength when Aq = 100 pm, but becomes 0.94 when 
Ao = 200 pm. Therefore, in order to fit the peak fiux density, the black body term will need to be smaller in the 
case of Ao = 100 pm, which means that the derived Tmbb must be smaller. For further demonstration, we generated 
a series of MBB models of varying Tmbb using Ao = 100 pm, convolved them with the SPIRE band response curves, 
and th en fi tted the simulated photometry using the MBB models with Ao = 200 pm. A representative case is given in 
Figure EH where the simulated objects are at 2 ) = 2 and all have = 10^^ L©. The left panel shows the relation 
of the derived Tmbb values and the input values. The peak temperature Tpeak^ on the other hand, will not be affected 
significantly because the fitting procedure, regardless of the choice of Ao, will always find the model that best matches 
the given SED, and therefore the peak wavelength and the peak temperature based on the Wien’s displacement law, 
will not change much. 

Similarly, the choice of Ao has little impact to the derived values because the shape of the best-fit model 

is governed by the observed FIR SED. When its peak region is well sampled by the three SPIRE bands, the small 
differences in the fitted SED beyond the peak will only have a small contribution to and thus have little influence. 
The right panel of Figure EH shows the comparison of based on the same simulation. 


B. SUMMARY OF ONLINE DATA TABLE AND EXAMPLES OF SED FITTING 

The header information of the online data table is given in Table IBll We note that this table includes all the 354 

sources as summarized in § 12.3.31 and Table [H _ 

The SEDs and the best-fit models for the IR quasars in the SNR5 subsample are shown in Figure EH as examples. 
For clarity, we only plot the MBB (red curves) and the SK07 (grey curves) model fits. While there are 134 objects 
in this subsample, only 102 of them have photometry in all the three SPIRE bands to allow for the MBB fit, and 
hence we only show these 102 objects in this figure. The other 32 objects in the SNR5 subsample still have Ljr values 
as determined using the starburst models (see which are included in the online data table. We note that seven 
objects in this figure also have PACS data, for which we also have MBB+PL fits (shown as the blue curves). 


C. TEST OF AGN/STARBURST DEGOMPOSITION 

The availabilit y of mid-to-far IR data in the rece n t years have allowed better characterization of AGN SEDs in 
this regime (e. g. iShi et all 120071: iNetzer et al.l 120071: iMaiolino et al.l 120071: iShang et al.l 120111: iMnllanev et al.l 120111: 
HhT et al.l [ 201 ^ iDale et all l2014f ). In light of these improvements, it has been advocated that the mid-to-far IR 
contributio ns from the AGN and the host galaxy star formation can be separated through the decomposition of the 
SEDs (e.g. IMnllanev et al.l [201 ll: iDale et ^l2014f ). Such decomposition approaches, however, still have caveats. An 
important one is that the AGN templates in use can only be constructed by subtracting the star formation contributions 
from the observed SEDs, where such contributions are determined through the calibration of some indicators of star 
formation activities (such as the PAH features). It is unclear whether such calibration is universally applicable given 
the complexity of the star forming regions. Furthermore, the far-IR part of the templates is just the extrapolation 
from the mid-IR part, and a cutoff has to be applied beyond a certain wavelength in FIR that is somewhat arbitrarily 
chosen. This means that such decomposition schemes are not appropriate to address the question of AGN contribution 
in FIR, because it is already assumed to be minimal by design. 

With these caveats in mind, we tested the AGN and the star formation (AGN/SF) decomposition of our objects in 

order to further check the consistency of our conclusion that is mainly due to the heating from star formation. We 
only carried out this test for the objects that have the PACS data and/or those that have the Spitzer MIPS 24/70 pm 
data readily available from the releases of the relevant Spitzer Legacy Survey programs residing in the IPAC Infrared 
Science Archinve (IRSA). In brief, there are in total 87 (COSMOS 33, Lockman 26, ELS 14, XMM-LSS-SWIRE 14) 
objects in our IR quasar sample that have at least MIPS 24 pm data, among which 37 objects have detections in all 
the three SPIRE bands, and thu s are suitable fo r the t est and the comparison. 

We used the software suite of IMnllanev et aP (j2Qlll ) in this test. This particular flavor of decomposition involves 
one AGN template and five starburst templates (denoted from “SBl” to “SB5”), and in each fitting run the routine 
combines the AGN template and one of the five starburst templates and determines the normalizations through the 
least-y^ fitting. The AGN template is forced to cut off at a wavelength between 20 and 100 pm in the form of a black 
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Table B1 

Summary of the online data table 


Column 

Name 

Description 

1 

SDSS_IAU_name 

SDSS quasar lAU name 

2 

HerscheLIAU_name 

lAU name of the Herschel counterpart as in the original catalogs 

3 

strJd 

Object ID assigned in this work 

4 

RA_250 

Herschel RA (J2000) 

5 

Dec_250 

Herschel Dec (J2000) 

6 

F250 

250 pm flux (mJy) 

7 

E250 

250 pm flux error (mJy) 

8 

F350 

350 pm flux (mJy) 

9 

E350 

350 pm flux error (mJy) 

10 

F500 

500 pm flux (mJy) 

11 

E500 

500 pm flux error (mJy) 

12 

is_SNR5 

1: SNR in 250 pm >= 5; 0: SNR in 250 pm < 5 

13 

RA.IOO 

100 iim RA (J2000) 

14 

Dec_100 

100 iim Dec (J2000) 

15 

FlOO 

100 pm flux (mJy) 

16 

ElOO 

100 pm flux error (mJy) 

17 

RA_160 

160 iim RA (J2000) 

18 

Dec_160 

160 iim Dec (J2000) 

19 

F160 

160 pm flux (mJy) 

20 

E160 

160 pm flux error (mJy) 

21 

RA_DR7 

SDSS DR7 RA (J2000) 

22 

Dec_DR7 

SDSS DR7 Dec (J2000) 

23 

z_DR7 

SDSS DR7 redshift 

24 

Mi_DR7 

SDSS DR7 z'-band absolute magnitude (at z=0) 

25 

RA_DR10 

SDSS DRIO RA (J2000) 

26 

Dec_DR10 

SDSS DRIO Dec (J2000) 

27 

z_DR10 

SDSS DRIO redshift 

28 

Mi_DR10 

SDSS DRIO z'-band absolute magnitude (at z=2) 

29 

LOG_MBH 

Black hole mass measured bv Shen et al. f2011j fLoe Mr;^j 

30 

LOG_MBHERR 

Error of the black hole mass 

31 

LOGXX 

Hard-band X-ray luminosity in the restframe 2-10 keV (erg/s) 

32 

GAMMA 

Photon index of the X-ray SED used to derive X-ray luminosity 

33 

XRAY_REF 

References of the X-rav data. CSC: Evans et al. f2010j: 3XMM: Rosen et al. f2015j 

34 

LOG_LIR_MBB 

measured using MBB fitting (Log Lq), integrated over 8-1000 pm 

35 

LOG_LIRERR_MBB 

Error of the measured 

36 

CHISQ_MBB 

of the MBB fitting 

37 

TMBB 

Derived black body temperature of the best fit model (K) 

38 

TMBBERR 

Error of Pmsu (K) 

39 

TPEAK 

Derived peak temperature of the best fit model (K) 

40 

LOG_MDUST 

Dust mass (Log Mq) 

41 

LOG_MDUSTERR 

Error of dust mass 

42 

LOG_SFR 

SFR usine: Kennicutt fl998L modified for Chabrier IMF fLoe: Mr:^/vri 

43 

LOG_MGAS 

Gas mass converted from dust mass (assuming gas-to-mass-ratio of 140) 

44 

LOG-TDEP 

Gas depletion time scale (Log yr) 

45 

LOG_LIR_SK07 

Li Pi measured using SK07 templates (Log Lq) 

46 

LOG_LIRERR_SK07 

Error of the SK07 Ljr 

47 

CHISQ_SK07 

of the SK07 fitting 

48 

LOG_LIR_CE01 

Lm measured using CEOl templates (Log Lq) 

49 

LOG_LIRERR_CE01 

Error of the CEOl Lm 

50 

CHISQ_CE01 

of the CEOl fitting 

51 

LOG_LIR_DH02 

Lm measured using DH02 templates (Log Lq) 

52 

LOG_LIRERR_DH02 

Error of the DH02 Lm 

53 

CHISQ_DH02 

X^ of the DH02 fitting 


body, and the exact cut-off wavelength is determined during the fitting. The exact AGN-hSB combination that gave 
the smallest was deemed to be the best fit, from which the total IR luminosities (integrating over 8-1000 pm) due 
to the AGN and the SB heating were then obtained. Here we denote these as and ^ respectively, 

and denote the fraction of the AGN contribution as fracf^^ = j -j- Figure ICT] compares 

thus obtained to (filled red circles) and (open blue circles) that we derived based on the fit to the 

SPIRE data as described in 0and 0 The mean difference between (—0.12dex) 

with a sample standard deviation 0.2 dex (0.2 dex). The results are fully consistent with what we derived using solely 
the SPIRE bands, i.e., and , which both measure the total IR luminosity due to star formation, agree 

to within ^ 17%, while L'f'^^ which measures only the total IR luminosity from the cold-dust component, is less than 
by 32%. Gonsidering the differences in the methods and the template sets employed, we conclude that the 

decomposition test also provides results that are consistent with our interpretation that is mainly due to the 

heating of star formation. 
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Wavelength (/im) 


Figure Bl. FIR SEDs and their best-fit models for the IR quasars in the SNR5 subsample. The dark blue points are the observed 
SEDs, while the gray and the red curves are for the best-fit SK07 and MBB models, respectively. For objects with PACS data, the best-fit 
MBB+PL models are shown in blue. 

D. STACKING ANALYSIS OF HERSCHEL-VNDETECTED SDSS QUASARS 

To investigate how the survey limits impact our results, we performed a stacking analysis of the SDSS quasars that 
are not matched in the adopted HerMES, H-ATLAS and HerS catalogs. While a fraction of these unmatched SDSS 
quasars could be blended sources rejected by our stringent matching criterion, the majority should be those that are 
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Figure Cl. Comparisons of derived IR luminosities using only the FIR SEDs to those obtained through the AGN/SF decomposition 

method of IMullanev et ahl (|2011h incorporating additional mid-IR data. For simplicity, only and are compared to 

here, plotted in red and blue, respectively. The left panel shows the differences with respect to the AGN fraction to Ljji (over 8-1000 pm) 

as determined by the decomposition, while the right panel shows these differences with respect to (effectively The red and 

the blue dashed lines show the mean values of the differences. The histogram of the differences are shown in the middle panel. 



Figure Dl. Stacked SPIRE image stamps of the Herschel-Mndetected quasars in the HerMES L6-XMM-LSS-SWIRE field at redshift 
2 < z < 3. The left panel shows the stamps generated using average in stacking, while the right panel shows those derived using median. 
In each panel, the stacks of the “invisible”, “weak”, “combined” and “contaminated” sets are shown from the top to the bottom. 

undetected in the current SPIRE 250 pm data. It is difficult to stack all such objects from all fields, because different 
fields have different survey sensitivities. Therefore, we only stacked on a field-by-field basis. We show here the results 
in the HerMES L6-XMM-LSS and L5-Bootes fields and the HerS field, which are representative of the entire data sets 
in terms of survey limit. 


D.l. Stacking procedure and photometry 

The unmatched SDSS quasars in each field were first split into different bins by redshifts. Eor simplicity, we only 
discuss the “equal redshift” case but not the “equal volume” one. We adopted a large step-size, Az = 1, to ensure 
enough number of objects for stacking at the most critical redshifts. The objects were then visually inspected in 
250 pm and were subsequently divided into three categories: ‘‘weak” objects seem to have some weak enhancement 
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Table D1 

Summary of photometry of stacked image 





invisible set 




weak set 



combined set 


z 


N 

‘S'250 

S 350 

‘S'500 

N 

S 250 

S 350 

‘S'500 

N 

S 250 

‘S'350 

S'soo 

Bootes-HerMES 

(0, 

1] 

24 

- 

- 

- 

10 

10.00 ±0.78 

2.53 ±0.78 

- 

34 

4.04 ± 0.54 

- 

2.48 ±0.64 



- 

2.48 ±0.68 

- 


9.32 ± 0.79 

- 

- 


4.73 ±0.54 

- 

- 

(1, 

2] 

21 

- 

- 

- 

15 

11.19 ±0.56 

7.17 ±0.56 

9.14 ±0.65 

36 

4.40 ± 0.47 

1.70 ±0.46 

1.68 ±0.51 


- 

- 

3.43 ±0.75 


11.80 ±0.56 

8.80 ±0.56 

7.48 ± 0.66 


4.76 ± 0.47 

2.53 ± 0.46 

3.23 ±0.55 

(2, 

3] 

66 

- 

- 

- 

22 

6.71 ±0.48 

5.59 ±0.48 

3.11 ±0.57 

88 

1.17 ±0.30 

- 

- 


- 

- 

- 


7.76 ± 0.48 

7.05 ± 0.48 

2.92 ±0.57 


1.67 ±0.30 

1.27 ±0.29 

- 

(3, 

4] 

21 

- 

- 

- 

4 

10.07 ± 1.18 

14.20 ± 1.19 

17.73 ±1.47 

25 

- 

- 

3.38 ±0.58 


- 

- 

- 


11.43 ± 1.19 

8.99 ± 1.16 

13.32 ± 1.52 


- 

- 

2.91 ±0.57 

(4, 

5] 

2 

- 

- 

18.48 ±2.00 

0 

- 

- 

- 

2 

- 

- 

18.48 ± 2.00 


- 

- 

18.48 ±2.00 


- 

- 

- 


- 

- 

18.48 ±2.00 

XMM-LSS-SWIRE 

(0, 

1] 

31 

- 

- 

- 

12 

9.44 ± 0.78 

6.37 ±0.79 

4.39 ±0.89 

43 

1.45 ±0.44 

- 

- 



- 

- 

- 


9.26 ± 0.78 

6.10 ±0.79 

3.85 ±0.89 


- 

- 

2.71 ±0.51 

(1. 

2] 

10 

- 

- 

- 

8 

12.47 ± 1.01 

6.85 ±0.99 

- 

18 

5.93 ± 0.74 

- 

- 


- 

- 

- 


12.45 ± 1.02 

9.00 ±0.99 

- 


6.47 ±0.74 

3.87 ±0.70 

- 

(2, 

3] 

113 

- 

- 

- 

41 

9.35 ±0.45 

5.92 ±0.44 

4.36 ±0.53 

154 

2.31 ±0.23 

1.70 ±0.22 

- 


- 

1.14 ±0.26 

- 


9.27 ±0.45 

5.45 ± 0.44 

4.76 ± 0.53 


2.71 ±0.23 

2.42 ± 0.22 

1.34 ±0.27 

(3, 

4] 

18 

3.28 ±0.83 

- 

- 

7 

9.16 ±1.12 

9.39 ± 1.05 

7.28 ±1.24 

25 

2.35 ±0.67 

- 

- 


- 

- 

- 


9.72 ± 1.12 

13.25 ± 1.05 

6.72 ±1.24 


3.63 ±0.67 

- 

2.98 ±0.72 

(4, 

5] 

1 

- 

- 

24.32 ± 2.93 

0 

- 

- 

- 

1 

- 

- 

- 


- 

- 

24.32 ± 2.93 


- 

- 

- 


- 

- 

- 

HerS 

(0, 

1] 

904 

7.12 ±0.63 

4.14 ±0.84 

2.22 ± 0.73 

0 

- 

- 

- 

904 

7.12 ±0.63 

4.14 ±0.84 

2.22 ±0.73 



6.88 ± 0.63 

4.26 ± 0.84 

2.53 ± 0.73 


- 

- 

- 


6.88 ±0.63 

4.26 ± 0.84 

2.53 ±0.73 

(1, 

2] 

699 

- 

- 

- 

452 

17.17 ±0.82 

11.22 ± 1.10 

6.67 ±0.99 

1151 

7.44 ± 0.54 

5.13 ±0.71 

2.82 ±0.63 


- 

- 

- 


17.35 ±0.82 

11.31 ± 1.10 

6.62 ±0.98 


7.66 ±0.54 

5.50 ±0.71 

3.07 ±0.63 

(2, 

3] 

1046 

- 

- 

- 

350 

14.97 ± 1.00 

11.54 ± 1.37 

8.46 ±1.15 

1396 

4.24 ± 0.53 

3.88 ± 0.70 

3.39 ±0.61 


- 

- 

- 


15.37 ± 1.00 

11.75 ± 1.37 

8.10 ± 1.15 


4.54 ±0.53 

3.84 ±0.70 

3.30 ±0.61 

(3, 

4] 

220 

- 

- 

- 

77 

13.11 ± 1.83 

10.46 ±2.52 

10.55 ±2.13 

297 

3.96 ±1.04 

4.79 ±1.41 

4.53 ± 1.27 


- 

- 

- 


13.31 ± 1.83 

10.91 ±2.52 

12.15 ±2.13 


3.99 ±1.04 

4.78 ± 1.41 

- 

(4, 

5] 

27 

- 

- 

- 

0 

- 

- 

- 

27 

- 

- 

- 
















Note. — Two algorithms, namely, average and median, were employed while stacking. The first row for each redshift lists the photometry 
using the average stacks, while the second row lists those using the median stacks. The fluxes are given in mJy; None detections are denoted 
using 


at the source locations, ‘‘invisihle^^ objects do not have any sign of detection, and “contaminated^^ objects are likely 
contaminated by nearby sources. The results from this classification , wh ile done in 250 pm, were assumed to be valid 
in 350 and 500 pm as well. The detailed numbers are listed in Table lun for these three fields. The stacking was done 
for these three categories separately in each redshift bin, using both average and median as the combining methods. 
We also combined the “invisible” and the “weak” objects into the “combined” set and stacked them together. In 
most cases, the “invisible” stacks still do not result in significant detections, while the “weak” stacks usually result 
in positive detections in at least 250 and 350 pm. The “combined” stacks, on the other hand, are weaker than the 
“weak” stacks, but are still positive detections in at least 250 pm. The “contaminated” stacks, as expected, often 
show contamination from the residuals of the blended neighbors, and hence were excluded from further discussion. 
For demonstration. Figure [DT] shows the average and the median 250 pm stacks in these four cases for the quasars at 
2 < ^ < 3 in the L6-XMM-LSS-SWIRE field. 

We used the Herschel Interactive Processing Environment (HIPE) software package to do photometry on the stacks. 
The Herschel science and uncertainty frames of the image cutouts at the undetected SDSS quasar coordinates were 
stacked separately and inputted to HIPE task sourceExtractorSussextractor, which then performed PSE fitting 
to measure the fluxes and asso ciat ed errors. A 3 a detection limit was adopted during the photometry. The measured 
results are also listed in Table [PTl 

As the “invisible” set usually does not result in detections, it will not be discussed further as an independent set. 
We will focus on the “weak” and the “comb” sets. 
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Table D2 

Summary of stacking analysis 




weak set 



combined set 



0 < 2 ; < 1 

1 < 2 ; < 2 

2 < ^ < 3 

3 < ^ < 4 

0 < z <1 

1 < 2 ; < 2 

2 < z <3 

3 < ^ < 4 

Bootes-HerMES 

^ 0 ^ 

T CEOl 

TBho2 

10 

15 

22 

4 

34 

36 

88 

25 

11.5 (11.5*) 
11.4 (11.4*) 

12.0 (12.0) 
12.0 (12.2) 

12.2 (12.3) 
12.2 (12.2) 

12.8 (12.7) 

12.8 42.6) 

11.2 (11.2*) 
11.2 (11.2*) 

11.5 (11.6) 
11.4 5i-7) 

11.5*(11.5) 

11.4*4i.7) 

12.2*(12.2*) 

12.2*(12.2*) 

XMM-LSS-SWIRE 

^ 0 ^ 

T CEfol 

TBh02 

^IR 

12 

8 

41 

7 

43 

18 

154 

25 

11.6 (11.5) 
11.4 51-4) 

12.1 (12.1) 
12.2 (12.2) 

12.3 (12.3) 

12.4 ^2.2) 

12.6 (12.7) 
12.6 ^2.6) 

10.5*(11.6*) 

10.6*(11.2*) 

11.8*(11.8) 

11.7*4i-7) 

11.7 (11.8) 
11.7 5i-7) 

12.1*(12.2) 

12.0*(12.2) 

HerS 

T CEfOl 

rBh02 

^IR 


452 

350 

77 


1151 

1396 

297 


12.2 (12.2) 
12.2 (12.2) 

12.6 (12.6) 
12.6 42.6) 

12.7 (12.8) 

12.8 ^2.8) 


11.8 (11.9) 
12.0 ^2.0) 

12.0 (12.0) 
12.0 ^2.0) 

12.3 (12.3) 

12.4 (12.2) 


Note. — The listed IR luminosities (LogL^^/LQ) with and without parenthesis are derived from the median and the average stacks, 
respectively. The asterisks denote that only one photometry point is available when doing the SED fitting. 

D.2. IR luminosity 

The FIR SEDs of the stacks could be fit in the usual way. There is a caveat, however. The implicit assumption 
of any stacking analysis is always that the sources being stacked can be scaled linearly. This is not necessary true in 
our case, because these Herschel-nndetected SDSS quasars could be of different dust temperatures and hence their 
intrinsic FIR SEDs could vary significantly in shape. Therefore, the physical properties inferred from the SED fitting 
of the stacks should be interpreted with caution. Here we only focus on the total IR luminosities of the stacks, and 
assume that these are the representative values of the stacked population. 

In most cases, only the 250 pm stacks have detections of sufficient S/N, and this prevented us from doing the full 
suite of FIR SED fitting as in ^ First of all, we had to abandon the MBB fit because it would require at least three 
bands to obtain reasonable constraints. This left the three starburst models to use, which all allow the estimation of 
^IR searching for the template whose associated Lir is the closest match to the one or two-band flux densities at 
the given redshift. However, we also opted to abandon the SK07 models because in some cases the results had large 
discrepancies with respect to those based on CEOl and DH02. This is largely due to the wide dynamic range and 
the fine grid of the SK07 templates, which could make their results highly degenerated in case when only one or two 
bands are available. As we only have a very limited number of bins here, including these deviant SK07 results would 
significantly skew the picture. Therefor e, w e used the results from the CEOl and the DH02 templates for this analysis. 
These results are summarized in Table [D2l 

D.3. Correction to fraction of IR-luminous quasars 

One key objective of this stacking analysis was to study how the distribution of the fraction of IR-luminous optical 
quasars (see Figure [15]) will change if the undetected population is taken into account. In particular, we are interested 
in understanding whether the sharp drop from z ~ 2 to higher redshifts is caused by the potential bias when using 
only the SPIRE-detected objects for the statistics. 

This IR-luminous fraction depends on the exact choice of Lc, which could be below the survey limit of the data, and 
therefore ideally we should co mpe nsate for the number of objects that are below the survey limit and yet still have 
^IR ^ ^c- However, as Table |D2] shows, this is not always possible because we cannot control the input undetected 
objects so that they will stack to the exact Lc required. Therefore, we need to make corrections in a diffe rent way. 

As we have chosen Lc = Lq, it is only necessary to analyze the “weak” sets, because from Table ID^ it is obvious 

that the individually undetected objects that still have Lj^ > Lc can only be found in these sets. We first start from 
the 2 < z < 3 bin. In the L5-Bootes field, there are 18 SPIRE-detected quasars in this bin that have Lj^ > Lq. 

The stack of 22 undetected objects in this bin, on the other hand, has the median of from the 

fits using CEOl or DH02 models. We use the superscript to denoted that this quantity is for the “undetected” 
objects. By the definition of median, half (i.e., 11) of these undetected objects should have Lj^ > Lq. As 

Lc > Lq, we conclude that at most there are 11 objects that can be added to the existing 18 objects. In other 

words, the number of objects above the threshold Lc after this correction, Ncor^ should be 18 < Ncor < 29, and 
the correction factor, /cor, should be fcor < 1-b. Similarly, for the bin of 1 < z < 2, there are 5 objects that have 
^IR ^ Lq, and the stack of 15 undetected objects has the median of Lj^'^ = 10^^’^ Lq. Therefore we get 

5 < Ncor < 12, and fcor < 2.4. 

We can repeat the similar analysis for all the fields, and obtain the constrained corrections for each redshift bin. If 
< I/c, such a correction is “maximum possible”, and if > Lc, such a correction is “minimum necessary”. 

At 1 < z < 2, the corrections are all “maximum possible”. At 2 < 2 ; < 3, the correction is “minimum necessary” for 
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the HerS field, and are “maximum possible” for all other fields. However, the HerS field dominates the statistics here 
at 2 < 2 ; < 3: the stack of 350 iJer5c/ie/-undetected quasars has median Lq, which means that we need 

to add 175 objects for the correction by this field alone. The total number of Herschel-detected quasars that have 
^IR ^ L 0 in this redshift bin is 84, and if we ignore the “maximum possible” corrections in other fields, the 

“minimum necessary” corrections will be at least (175 + 84)/84 ^ 3.1 for the entire sample. 

D.4. Correction to IR luminosity density 

As expected, the Herschel-nndetected SDSS quasars are significantly less luminous than the Herschel-detected 
ones at the same redshifts. However, their collective contributions to are comparable to those from the 

Herschel-detected population, simply because of their larger number. This is particularly true at z = 2-3, where their 
contribution is even higher by a factor of ^ 2 x. 

The correction in each redshift bin was done by adding the product of the average stacked Lj^ values of the 
“combined” sets multiplied by the number of input objects in these sets. After such corrections, the evolution of 
will peak at z ^ 2.5 instead of z ^ 1.5. 
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